高速sqrt()関数


皆さんもご存知でしょうが<libps>標準の<sqrt()関数>は遅い.
ふざけんなってぐらい遅い.
解析しようかと思ったんですけど,面倒くさいし長ったらしいのでやめました.
まあ,大きな要因はPlayStationが浮動小数点の演算をソフトウェアでやっていることが原因なんでしょうけど.
PlayStationで浮動小数点はあまり使わないほうがいいようです.
なぜコプロセッサをつまなかったのだろう?
そんなに使うものでもないだろうけど....

まずは平方根の計算のアルゴリズムから考えよう.
早速公式から見ると,

  X
n+1 = (Xn + X1 / Xn) / 2
  (1 < n , n = ∞)

この公式で近似値を計算することができます.
実際に<2>の平方根を求めてみると,

  X
1 = 2
  X
2 = (X1 + X1 / X1) / 2 = (2.0000 + 2 / 2.0000) / 2 = 1.5000
  X
3 = (X2 + X1 / X2) / 2 = (1.5000 + 2 / 1.5000) / 2 1.4167
  X
4 = (X3 + X1 / X3) / 2 = (1.4167 + 2 / 1.4167) / 2 1.414215744335427

本当の答えは「
ひとよひとよにひとみごろ」で,1.414213562373095….

  近似値 : 1.41421
5744335427
  正解値 : 1.41421
3562373095

なんとほとんど一緒じゃねーか!!
3での計算結果の小数部の桁をもっと増やせばさらに近くなるはず.
または,X
5・・・Xnと計算するともっと近くなる.

ここで1つ謎があるんだけどね,<n>はどうやって決定するのか?
nが大きければ大きいほど,近似値は良い精度になるわけです.
でも,永遠にやっていたらきりがない上に標準のものよりも遅くなっちまう恐れがある.
そこで,小数点の桁を決めておき,次回の計算も同じになったら終了させる.

  X
1 = 2
  X
2 = (X1 + X1 / X1) / 2 = (2.0000 + 2 / 2.0000) / 2 = 1.5000
  X
3 = (X2 + X1 / X2) / 2 = (1.5000 + 2 / 1.5000) / 2 1.4167
  X
4 = (X3 + X1 / X3) / 2 = (1.4167 + 2 / 1.4167) / 2 1.4142
  X
5 = (X4 + X1 / X4) / 2 = (1.4142 + 2 / 1.4142) / 2 1.4142

上の計算でX
4 = X5なので,ここで終了させる.

こんどは高速化を考えてみる.
PlayStationでは浮動小数点はサポートされているものの,
ほとんどの計算は固定小数点で行われているように思える.(ほんとかな〜?)
それに浮動小数点での演算は処理が遅いので固定小数点がいいに決まってる.
これで,標準よりもかなり早くなる.
以下,サンプルソース.

  u_long sqrt2(u_long f)
  {
    u_long s = f,t;
    
    if(x == 0) return 0;
    do
    {
      t = s;
      s = (t + f / t) >> 1;
    }while(s < t);
    return t;
  }

  ※引数の小数部bit数 >> 1 = 戻り値の小数部bit数
   例えば,引数の小数部を24bitとした場合,戻り値の小数部は半分の12bitになります.

しかし,問題があるんですよ.
それは,精度が少し落ちてしますことです.
例えば,<2>の平方根を求めた場合

  sqrt(2.0)   = 1.414
214
  sqrt2(2 << 24) = 1.414
063
  ※固定小数点の小数部を24bitとする

この場合の処理時間を計測すると,

  sqrt(2.0)   ⇒
136
  sqrt2(2 << 24) ⇒
1
  ※ルートカウンタを使用して計測

精度をもっと良くしたい場合,以下のようにソースを変更する.

  float sqrt2(float f)
  {
    float s = f,t;
    
    if(x == 0) return 0;
    do
    {
      t = s;
      s = (t + f / t) / 2;
    }while(s < t);
    return t;
  }


この場合,
標準のものと結果は同じになります.(多少違いはあるかも!?)
ですが,浮動小数点を使っているのでやっぱり処理速度は落ちます.

  sqrt(2.0) ⇒
136
  sqrt2(2.0) ⇒
89
  ※ルートカウンタを使用して計測


結果的に,
どちらをとるかというところでしょうな.
  1.ダントツな速さだが,精度が落ちる方法.
  2.精度はいいが,処理速度が落ちる方法.

どちらも処理速度を上げようと思ったら,あとはアセンブラでのコーディングぐらいかな?
<1.>の精度を上げるには小数点部を多くとると良いでしょう.
私としては,<1.>の方法がお勧めですね.
処理速度がけた違いに速いし,精度もそれほど気にするものでもないように思えるから.
ちなみに,<1.>をアセンブラでコーディングしてみたので,良かったら使ってみてください.
こちらからどうぞ
処理速度はCでコーディングしたものとそんなに変わらないんですけどね.
以下,計測結果.

  C     ⇒ 145
  アセンブラ ⇒
126
  ※<2>の平方根の計算を1000回実行した結果
   ルートカウンタを使用して計測

見よ!!
1000回行ってもこの速さ.シビレるね〜.

Back