tsuy and buyobuyon - ttth @ ウィキ

コンパイラ組み込み関数

tsuyさんのページにいろいろ書き込んでいるうちに書き込める量の限界に達してしまったので、
コンパイラに組み込んだ浮動小数点数演算関数に関して、以前私が書き込んだ部分を
こちらに移動しました。

• ええと、ライブラリの件なのですが、
  まず、浮動小数点数演算の
  
  atan : type(fun type(float) -> type(float))
  cos : type(fun type(float) -> type(float))
  sin : type(fun type(float) -> type(float))
  
  は完全にお任せします(とりあえずはもうできているのでしたよね)。
  それから、
  
  int_of_float : type(fun type(float) -> type(int))
  float_of_int : type(fun type(int) -> type(float))
  fless : type(fun type(float) * type(float) -> type(bool))
  
  は、FPUで実装されているので、コンパイラの拡張ですぐ対応できます。
  それから、
  
  fiszero : type(fun type(float) -> type(bool))
  fispos : type(fun type(float) -> type(bool))
  fisneg : type(fun type(float) -> type(bool))
  fabs : type(fun type(float) -> type(float))
  fneg : type(fun type(float) -> type(float))
  fhalf : type(fun type(float) -> type(float))
  fsqr : type(fun type(float) -> type(float))
  floor : type(fun type(float) -> type(float))
  
  は、簡潔なので、私が入出力命令と一緒にアセンブリ言語で記述するなり、
  コンパイラ内部での変換なりで対処しておきます。
  それらに関して、ちょっと仕様の確認をしたいのですが、
  fisposとfisnegは0レジスタとの比較で良いですよね。
  fabs/fnegは基本的には最上位ビット0/反転で良いと思うのですが、
  引数が0のときはどうしましょう?
  今のFPUでは、いわゆる-0が入力された場合に結果はどうなります?
  あと、fiszeroはfabs→0レジスタとの比較で、
  fsqrは同じ値を引数としたfmulで、
  floorはftoiとitofの合成で良いですよね。
  fhalfは…、÷2.0にするのは遅すぎる…?
  仮数部のみデクリメントってもっと厄介ですよね…。
  
  あと、浮動小数点数演算以外に関しては、
  
  xor : type(fun type(bool) * type(bool) -> type(bool))
  read_float : type(fun type(unit) -> type(float))
  read_int : type(fun type(unit) -> type(int))
  print_int : type(fun type(int) -> type(unit))
  
  いずれも対応する命令が存在するので、コンパイラの拡張ですぐ対応できます。あと、
  
  print_char : type(fun type(int) -> type(unit))
  
  が何者だかいまいち把握していないのですが…。
  wordアラインメントだからprint_intと一緒にするしかありませんよね…。 -- buyobuyon (2007-03-02 15:49:50)
• fnegは0から引けばよかったのでしたね…。
  コンパイラのアセンブリコード生成部分を作成したときにはそういう変換をしていたのに…。
  
  absは、もし最上位ビット反転という方針がだめな場合、
  例えば、fnegしたものと、もとの値のandをとるという方針なら大丈夫ですよね? -- buyobuyon (2007-03-02 15:56:36)
• もうひとつボケをかましていましたね…。
  fhalfは÷2.0より×0.5ですね…。
  
  あともう１つ。現在、コンパイラでは、
  fequalを通常の整数の比較で置き換えてしまっているのですが、
  これは大丈夫ですよね? -- buyobuyon (2007-03-02 16:09:26)
• >fabs/fnegは基本的には最上位ビット0/反転で良いと思うのですが、
  >引数が0のときはどうしましょう?
  >今のFPUでは、いわゆる-0が入力された場合に結果はどうなります?
  FMULでは+0が出力されます。
  FADD,FSUBでは+0と-0の区別はつけず、また、出力が0となる場合は必ず
  +0になります。（多分…）
  FLESSはご存知の通り-0<+0の判定に使います。
  FTOIは+0の場合と処理は同じです。
  FINVSQRTは非正規化数の処理がいらないということなので、挙動は未定義です。
  FPUの出力の0はすべて+0となっているはずなので、FLESS絡みで問題がなければ
  fabsは最上位ビット0でいいような気がします。
  fisneg(-0.0)はtrueでいいんですかね？ -- tsuy (2007-03-02 18:35:21)
• >floorはftoiとitofの合成で良いですよね。
  floorって確か小数点未満切捨てですよね？
  FTOIの結果は小数点未満四捨五入としているので問題になる気がするのですが、
  FPUのFTOIを切り捨てにしたほうがいいですかね？
  (そっちのほうが速くなるので好都合だったりしますが)
  
  >print_char : type(fun type(int) -> type(unit))
  >が何者だかいまいち把握していないのですが…。
  すみません、僕もよく分かりません…。 -- tsuy (2007-03-02 18:42:40)
• >absは、もし最上位ビット反転という方針がだめな場合、
  >例えば、fnegしたものと、もとの値のandをとるという方針なら大丈夫ですよね?
  なるほど、それなら大丈夫ですね。
  
  >fequalを通常の整数の比較で置き換えてしまっているのですが、
  >これは大丈夫ですよね?
  大丈夫だと思います。 -- tsuy (2007-03-02 18:54:22)
• fequalの話で思いついたんですが、丸めを切捨てにしてしまったので、
  (1.0/3.0)*3.0と1.0が等しくならないんですよね…
  FMULだけでも丸めを切り上げにしたほうがいいだろうか… -- tsuy (2007-03-02 18:58:12)
• なるほど…。+0.0と-0.0がきわどいですが、とりあえず、
  
  fiszero → rd := (rs AND 0x7fffffff) >. +0.0
  fispos → rd := rs >. +0.0
  fisneg → rd := rs <. -0.0
  fabs → rd := rs AND 0x7fffffff
  fneg → rd := +0.0 -. rs
  fhalf → rd := rs *. +0.5
  fsqr → rd := rs *. rs
  
  としてみます。
  floorをどうするか…。
  floorって切り捨てというよりも、～を超えない最大の整数を表すfloat値を
  返す演算でしたよね、正負で非対称な…。
  コンテストルールがあるので、ftoiの仕様を変えるわけにもいきませんし…。一応、
  
  floor → r1 := ftoi(rs); r1 := itof(r1); rd := rs <. r1; rd := itof(rd); rd := r1 +. rd
  
  とすると、正しい値が求まりそうですが(そうでもない？)、命令数が…。
  あとは…、0.5引いて四捨五入、みたいなアバウトなことが許されるのであれば
  速くなりそうですが…。何かいい方法はないですかね…。
  
  あと、FPUで提供する演算、やっぱりFINVSQRTではなくてFSQRTの方にしてもらえます？
  たぶんその方が高速かつ楽なので(特にFINVSQRTをFSQRTとFDIVの合成として求めているのであればなおさら)。
  FINVSQRTのライブラリでの実装も任意で良いようです。
  何か、たくさん実装してもらった後で申し訳ないのですが…。 -- buyobuyon (2007-03-02 20:30:25)
• FSQRTですか…
  以前開平法で実装したものはサイズ、速度ともに実用的ではなかったので、
  また一から作ることになると思いますが、コンテストまでに間に合うだろうか…
  それからサイズの方も、今のFINVSQRTを超えると完全にオーバーしてしまう
  という不安があります。
  
  しかしどうやって実装しよう… -- tsuy (2007-03-02 22:08:00)
• >FINVSQRTをFSQRTとFDIVの合成として求めているのであればなおさら
  すみません、これはどういうことでしょうか？ -- tsuy (2007-03-02 22:12:34)
• FPUにFSQRTをのせる場合、逆数や割り算はどうなります？
  ライブラリで作ったINVでは少々遅いので、やはりFPUには
  FINVSQRTがいいような気がしないでもないのですが… -- tsuy (2007-03-03 15:56:18)
• >>FINVSQRTをFSQRTとFDIVの合成として求めているのであればなおさら
  >すみません、これはどういうことでしょうか？
  あれ？私、何か勘違いしていました？
  ↓に
  >・FINVSQRT
  >　－SQRTとFDIVで計算
  と書いてあったので、
  FINVSQRTは、内部でFSQRTとFDIVを組み合わせているだけなのかなと思ったのですが。
  だとすれば、FSQRTの方がFINVSQRTよりは速いわけですよね。
  
  min-rtでは、SQRTのみ使用されていて、INVとINVSQRTは直接使用されてはいないので、
  SQRTではなくINVSQRTがFPUで提供された場合には、
  コンパイラ内部では、"1 ./ SQRT(x)"をINVSQRT(x)に変換し、
  それ以外の"SQRT(x)"をINV(INVSQRT(x))に変換することで高速化を図るわけですが、
  よくよくmin-rtを見ると、
  "1 ./ SQRT(x)"を求める場合でも、直前or直後に別の用途で
  "SQRT(x)"を使用している場合が多いのですよ。
  例えば、min-rtには、
  
  let l = sqrt (fsqr v.(0) +. fsqr v.(1) +. fsqr v.(2)) in
  let il = if fiszero l then 1.0 else if inv then -1.0 /. l else 1.0 /. l in …
  
  という部分があって、まぁfiszero lは、fiszero (fsqr v.(0) +. fsqr v.(1) +. fsqr v.(2))に変えても
  一般的に問題はないはずですが、一応意味の違う文面ですし、またコンテストルールに
  
  >非最適化実行の場合と同一オペランドの浮動小数点演算が
  >同一回数行われなければならない。ただし、(2)については定数
  >割り算最適化(後述)、冗長性除去、定数畳み込みは例外とする
  
  とあるので、この変換はちょっときわどいかなと思いまして。
  
  とすると、結局INVSQRTしたものだけではなくSQRTしたものも求める必要があります。
  
  INVSQRT(x)とSQRT(x)の両方を求める必要があるとしたら
  SQRTがFPUで求まったいた方が速く(速いですよね？)、
  またプログラムの文面そのままで変換すれば良いので私の方も楽にはなるのですが…。
  
  あ、でも、上の話はFPUでFINVSQRTは内部でFSQRTとFDIVを組み合わせることで実装している、というのを前提にしていたので、
  そうでないのであれば実装し直すのも大変ですし、このままFINVSQRTでやってみます。 -- buyobuyon (2007-03-03 17:14:42)
• なるほど、分かりました。
  
  下のほうでFSQRTとFDIVの合成となっているのは、11月にそのようにして実装したものが
  実用的な性能ではなかったのでニュートン法で実装しなおしたという経緯があるのですが、
  その際にwikiの記述を直し忘れていたという私のミスです。すみませんでした。
  
  FSQRTは開平法では前述の通り性能が悪く、ニュートン法では割り算をする必要があるため
  どうしても遅くなるので、ニュートン法のFINVSQRTを入れたほうが良いよねというお話を
  した記憶があります。また、FSQRTを採用するとFPUのサイズ的にINVをライブラリで実現
  せざるをえないと思うのですが、それだとやはり遅くなる気がします。FINVSQRTを使えば
  その結果を2つ掛けるだけですし。
  
  時間的余裕があればFSQRTもやってみたいのですが2週間では多分無理なので、FINVSQRT
  を使う方向でお願いします。 -- tsuy (2007-03-03 19:04:26)
• あ、了解です。
  私が勘違いしていただけなので、FSQRTは実装し直さなくて良いですよ。
  
  なら、結論として、
  FPUで実装されるのは、FADD/FSUB/FMUL/FINV/FINVSQRT/FTOI/ITOF/FLESSの8つで、
  コンパイラの方では、
  
  (x ./ y)　→　(x *. (inv y))
  (sqrt x)　→　(inv (invsqrt x))
  
  とし、(実装する時間があれば)最適化のステップで
  
  (inv (inv x))　→　x
  
  という変換をすれば良いのですよね。
  
  コンパイラは直接の関係はないですが、OPCODEってどうなってます？ -- buyobuyon (2007-03-03 19:36:24)
• FPUで実装されるのはFADD/FSUB/FMUL/FINVSQRT/FTOI/ITOF/FLESSの7つです。
  上のほうで誤解させてしまう書き方をしてしまいましたね。すみません。
  
  FINVは、すぐ下にあるとおり基盤の18×18の乗算器の数が足りないというのと、
  他の命令のパイプライン化でサイズが大きくなることが予想されるため、FPUに
  載せることは不可能です。
  
  だからおそらく
  (x ./ y)　→　(x *. （(invsqrt y) *. (invsqrt y)))
  (sqrt x)　→　(invsqrt （(invsqrt x) *. (invsqrt x)))
  のような感じになるかと思われます。
  
  OPCODEはFPU外部から下の仕様にあるFUNCの部分に4bit入力してもらうように
  しています。命令セットのページにあるOPCODEでFPUでサポートしている命令の
  OPCODEなら結果を規定のクロック後に出力し、それ以外なら
  "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
  が出力されるようになっています。 -- tsuy (2007-03-03 21:01:24)
• あ、了解です。
  割り算にもinvsqrtが入りますか…、結構遅延が大きくなりそうですね。
  FINVが入らないってことは、4bitずつ引き算っていうFDIVも入れるのは無理なのですよね？
  
  あと、FINVSQRTを使ってFDIVやSQRTを実現する場合、
  FDIVやSQRTの精度の条件は満たされますか？
  
  何かいろいろと聞いてしまって申し訳ありません。
  
  あと、FPUのOPCODEの件、一応その4bit値がわかるように、
  命令セットのページにでもまとめておいておいてもらえます？
  アセンブラとシミュレータも多少変更する必要があるでしょうし。 -- buyobuyon (2007-03-03 21:54:35)
• FDIVどころかFINVSQRT自体入る保証がなかったり…
  パイプライン化しなくていいのなら余裕なんですけどね。
  ちなみにFDIVはFINVSQRT10個分のサイズなので、使い物になりません。
  
  精度は…満たされないような気がしてきました。
  FDIVとSQRTに使うFMULとFINVSQRTの誤差は最大で仮数部23bit目のずれ
  なのはFPUのテストで分かっていて、FDIVもSQRTも4回これらが使われて
  いるんですよね…少し考えて見ます
  
  OPCODEの件、命令セットのところに書かれているOPCODEをそのまま4bit
  値にしているんですがそれをまとめるということですか？ -- tsuy (2007-03-03 23:12:08)
• 精度の件について
  FDIVですが、
  (x ./ y)　→　(x *. （(invsqrt y) *. (invsqrt y)))
  を使うとすると、
  (invsqrt　y)＝(y^(-1/2))*(1+2^(-23))
  であり、
  (invsqrt y) *. (invsqrt y)＝(invsqrt y)*(invsqrt y)*(1+2^(-23))
  なので
  (invsqrt y) *. (invsqrt y)＝(y^(-1/2))*(y^(-1/2))*（(1+2^(-23))^3)
  となりますよね？
  さらにこれにxを掛けるので、真の値との誤差は
  (1+2^(-23))^4　-　1
  分になるものと思われます。
  そしてこの値を計算すると
  1+2^(-21)+2^(-44)+2^(-45)+2^(-67)+2^(-92)
  となります。
  ルールは確か2^(-20)だったと思いますが、これならぎりぎりセーフかと思います。
  計算が違ったらすみません -- tsuy (2007-03-03 23:56:13)
• SQRTについてですが、
  (sqrt x)　→　(invsqrt （(invsqrt x) *. (invsqrt x)))
  を使うものとして、
  (invsqrt x) *. (invsqrt x)
  部分の真の値をaとおきます。
  上と同様に考えると、誤差を入れた値は
  a*（(1+2^(-23))^3)
  となり、これにFINVSQRTを適用すると、誤差も入れて
  [｛a*（(1+2^(-23))^3)｝^(-1/2)]*(1+2^(-23))
  =｛a^(-1/2)｝*｛(1+2^(-23))^(-1/2)｝
  となります。ルールでは確か2^(-20)だったかと思うので、
  |1　-　(1+2^(-23))^(-1/2)| < 2^(-20)
  なら問題ないことになります。
  実際計算してみると、同値変形で成立した不等式になったので、恐らく
  問題ないかと思います。
  
  とは言っても計算に自信がないのですが、他の高速化した命令が正しいか
  テストするときに一緒にテストしたほうがいいかな？
  ついでに、というには時間と手間がかかるけれど。 -- tsuy (2007-03-04 00:19:38)
• 大丈夫そうですね、それなら先ほどの式を使わせていただきます。
  
  ただ一つ気になったのが、FDIVやSQRTを求める際に(INV x)を求める段階で、
  let t = (invsqrt x) in (t *. t)
  とすると最初の演算の誤差が後の演算の2つのオペランドに影響しますが、
  let t = (x *. x) in (invsqrt t)
  とすると後の演算のオペランドが1つなので、
  結果として浮動小数点数演算1つ分誤差が小さくなるような気がするのですが、どうですかね？
  sqrtの入力が負の場合が未定義で良いのなら、ルール上は問題ないはずですし。
  (また何か勘違いしているような気もしますが…)
  
  テストは余裕があればで良いと思いますよ。
  結果的にレイトレが動けば良いわけですし(早く一度動かさないとね…)。
  
  あと、OPCODEの件、またまた私の勘違いでした。すみません。
  FINVが無いなら、前に定義したもので足りるのでしたね…。 -- buyobuyon (2007-03-04 00:29:59)
• 話が変わるのですが、floorの実装をライブラリの方でやっていただけますか？
  よくよく考えると、ftoiしてからitofする方針だと、intの範囲を超える整数のfloat表現が
  入力に入ったときに厄介ですし、
  まじめにやるとなるとマクロというより関数という感じのサイズになってしまうもので。
  一応、関数内で何をすればよいかはわかっている(つもりな)ので、
  もし大変であれば私が実装しておきますが、どうしましょう？ -- buyobuyon (2007-03-04 00:56:14)
• 確かに
  let t = (x *. x) in (invsqrt t)
  の方が誤差が少なくなりますね。盲点（？）でした。
  
  そちらのやることもいっぱいあるようですし、floorの実装をさせていただきます。
  ただ性能のよくないものを作ってしまってはいけないので、方針を教えていただけますか？
  それからもう一つ、OCamlの文法で書いていいのかな？ -- tsuy (2007-03-04 01:54:20)
• すみません、お願いします。
  
  実装はOCaml文法で構いません。できれば、アップしたコンパイラで
  コンパイルできる形にすぐ直せるよう、多相関数などの使用は控えてください。
  あと、ftoiはint_of_floatという名前の一変数関数として、
  itofはfloat_of_intという名前の一変数関数として、それぞれ実装しますので、
  必要なら使ってください。
  
  ええと、実装方法についてあれから改めて考えてみたのですが、
  floatの仮数部のsizeってsizeof(int)より断然小さいですよね。
  そのため、floatで表された値がintで表現できる範囲を超えた場合というのは
  もう既にfloatの仮数部に1を足したら(引いたら)float値は1以上変化する場合なので、
  そのときにはそのまま引数を返してやればよいわけですよね。
  
  とすると、指数部がその境界となる値以上であったら引数をそのまま返し、
  そうでなければ先ほど述べたように、
  FTOI → ITOF → FLESS結果が真なら(int値に)-1 → ITOFというので
  良いかなと思うのですが、どうでしょう？
  たぶん、下手に指数部に応じてシフト・場合分けを繰り返すよりは
  このようにFTOI→ITOFでごまかした方が速いような気がするのですが…、
  あまりクロック数とか意識していなかったので、どうなんだろ？
  
  一応、上の方針だと2箇所だけ場合わけがありますが、2つ目に関しては、
  真=int型の-1なので足し算で代用可能ですよね。
  とは言っても、ML文法だとint値にbool値を足すようなことはできないので、
  後で自作コンパイラでコンパイルして出力されたコードを一部書き換えるように
  すれば良いと思います。
  他のcosやatanなどもその方針で高速化できるものがあるかもしれませんし。
  
  あと、finvは結局let t = (x *. x) in (invsqrt t)の形を採用しました。 -- buyobuyon (2007-03-04 03:15:41)
• なるほど分かりました。
  
  念のため確認しておくと、
  入力ｘについて、
  ｜ｘ｜≧2^(23)のとき、この浮動少数はすでに整数なのでそのまま返す
  こうすることでftoiでINT_MAXやINT_MINが返る問題を避ける
  ｜ｘ｜＜2^(23)のとき、y=itof(ftoi(x))とすると、yはxを四捨五入した
  ものになるので、x＜yなら切り上げられているのでy-1を、
  x≧yなら切り捨てられているのでyを返す
  ということでいいですよね？
  
  上の方で
  >0.5引いて四捨五入、みたいなアバウトなことが許されるのであれば
  >速くなりそうですが
  とありますが、できないですかねー…
  出来るとありがたいんですけどね -- tsuy (2007-03-04 16:24:53)
• itof・ftoiによる方針、それで良いと思います。
  
  0.5引いて四捨五入も実はOKのような気がしてきたのですが…。
  0.5が非常にきりの良い値なので、引き算したときに滅多に桁落ち(情報落ち？用語がわからない…)しませんから、
  0.5の引き算で正しい値が求まる下限と上限の2つの境界で場合分けすれば、この方針でもいけますよね？
  ただ、場合分けが1つ増えると命令は2つ増えるので、
  結局命令数は上の方針と変わらないような気がします。
  
  あと、-0.5を四捨五入したときに0.0でなく-1.0になるように、
  すなわち正負で対象になるように四捨五入が定義されていると、
  後者の方針では正しい値は求まりませんね…。 -- buyobuyon (2007-03-04 17:15:28)
• ftoiによる四捨五入は正負で対象となっているので、0.5引くのはよくない
  ようですね。
  
  最初の方針で実装しようかと思います。
  テストも含めて火曜日夜あたりまでにこのページに上げようかと思うのですが、
  それで間に合いますよね？ -- tsuy (2007-03-04 19:12:55)
• あ、大丈夫だと思います。
  というか、私が火曜日夜までWiki見られなくなりますので…。
  一応こちらもテスト作業だけは進めておきます(できればある程度の最適化も).
  -- buyobuyon (2007-03-04 21:01:06)