math/04 のバックアップ(No.2)

[ トップ | 一覧 | 単語検索 | 最終更新 | バックアップ | ヘルプ ]

• バックアップ一覧
• 差分 を表示
• 現在との差分 を表示
• ソース を表示
• math/04 へ行く。
  • 1 (2007-04-01 (日) 10:43:43)
  • 2 (2007-04-01 (日) 11:33:18)
  • 3 (2007-04-01 (日) 11:33:18)

---

簡単な開根アルゴリズム

• (by K, 2007.04.01)

(1) ルートを計算したい

• アセンブラで遊んでいると、整数演算だけでsqrtの計算をしたいと思うことがある。その場合、I.Tak.さんみたいに筆算での開根アルゴリズムをうまく実装するのが正しいのだけれど、そこまでまともなものいらないっていう場合がある。これはそんなときにアルゴリズム。
• 中学生に教えてあげるときも使えるかもしれない。だって簡単だから。

(2) 基本原理

• 1から順番に奇数を足していくと次々と平方数ができる。以下に例を示そう。
  • 1 = 1^2
  • 1 + 3 = 4 = 2^2
  • 1 + 3 + 5 = 9 = 3^2
  • 1 + 3 + 5 + 7 = 16 = 4^2
  • 1 + 3 + 5 + 7 + 9 = 25 = 5^2
  • これは一般には n^2 + (2n+1) = (n+1)^2 という式の性質を利用しているので必ず成り立つ。
• これを逆に考える。任意の正の値xに対して、奇数を順番に引いてみる。引けた回数がsqrt(x)の値である。
  • x=10の場合の例：
    • 10から1を引いてみる。引けた！答えは9。
    • その9から今度は3を引いてみる。引けた。答えは6。
    • その6から今度は5を引いてみる。引けた。答えは1。
    • その1から今度は7を引いてみる。引けない！
    • 引き算が成功した回数は3回。だからsqrt(10)=3。
• （このアルゴリズムは僕が考えたものではない。中学か高校のときに本で読んだものである。）

(3) 精度

• sqrt(10)の結果が3であるというのは、いかにも精度が悪いかもしれない。もうちょっと下の桁が知りたいこともある。そうであるなら、xを100倍すればよい。
  • sqrt(x * 100) = sqrt(x) * 10
• そうすると実際に求めたい値の10倍が得られる。つまり31が得られる。これを10で割れば3.1、つまりsqrt(10)=3.1と分かる。精度が1桁改善した。
• 100倍ではなく、1万倍、100万倍すれば精度はもっと良くなる。
  • アセンブラで実装するときは、100倍とかではなく、4、16、64、256倍などを用いるほうがいいと思われる。

(4) 速度

• これで精度の問題も解決した。原理上は近似が入っていないため、これでどこまでも精度を上げられる。アルゴリズムも簡単でこれは最高！・・・と思ったらそれは間違いである。最大の問題として、このアルゴリズムは演算速度が遅い。sqrt(3)を1.732まで計算しようとしたら、sqrt(3000000)を計算することになるわけだが、これの結果は当然1732であって、つまり1732回も引き算をさせられるということなのである。これはいかにも重い。
• これを改良する方法としては、最初に1から引くのではなく、最初に1万を引いてみる、という手がある。もし引けたら、次は201から引いていくし、引けた回数ももう100回分は完了したとみなす。これで最初の100回分の引き算を省略できる。
• もし1万が引けた場合、さらに3万を引いてみてそれで引けたら、次は401、200回完了、とする。・・・とまあ、こうやって工夫すればそれなりに速くもなる。

(5) 単純な2分検索との比較

• こういう工夫をやりだすと、単純な2分検索と比べてどうなのかが気になってくる。単純な2分検索というのはこういうものだ（ここではintが32bitのCPUで面倒な多倍長処理をしないことを想定している）。

  unsigned int isqrt(unsigned int x)
  {
      unsigned int a = 0, b = 65535, c;
      do {
          c = (a + b) / 2;
          if (x < c * c)
              b = c;
          else
              a = c;
      } while (a + 1 < b);
      return a;
  }

• つまり最初に0〜65535という値の範囲を考える。そしてその真ん中の数32767を取り出して自乗し、xと比較する。xのほうが小さいようなら、範囲を0〜32767に絞る。xのほうが大きいようなら範囲を32767〜65535に絞る。これを何度か繰り返すうちに範囲が狭くなり、十分に狭くなったらやめる。これなら、ループ一回につき値の範囲は半分になっていくから、答えが16bit未満ならたった16回程度のループで答えが出せそうである。
• （なんかこれで十分かもしれない気がしてきた・・・。）

(6) 改良

• これで終わってしまうとあまりに悲しいので（笑）、もうちょっと考えてみることにする。
• sqrt(3)について、もっとたくさんの桁を求められないだろうか。つまり17bit以上を。
• sqrt(x * 100) = sqrt(x) * 10 なのだから、こんな方法はどうだろうか。
• 最初に100倍しないで計算し、sqrt(3)=1という結果を得る。しかしこれは本当は、sqrt(3)が1以上2未満であるということを意味している（切り捨てなので）。ちなみにこのときの引き残りは2だ。
• そこで精度を上げるために、sqrtの中を100倍し、答えも10倍する。そうすると、さっきの結果だけで計算しなくてもsqrt(300)=10という意味に取れる。もっというと、引き残り200の状態で回数10回から再開できる。次に引く数は21だ。200から、21, 23, 25, 27, ...と引いていくと33まで引ける。つまり追加で引けた回数は7回だから、合計すると17回。だからsqrt(300)=17とわかる。引き残りは11だ。
• さらに精度を上げるために、11を100倍し、回数170回から再開してみる。1100から341, 343, 345, ... と引いていく。345は引けるが、引いた結果は71なので、ここでおわる。ここまでの計算でsqrt(30000)=173とわかる。
• もう分かってきたと思うがあと一回だけやろう。さらに精度を上げるために7100の引き残りに対し、回数1730回からはじめる。3461, 3463とここでもう次は引けない。sqrt(3000000)=1732。引き残りは176。
• さてとりあえずここで打ち切るとすると、sqrt(300万)=1732という結果を出すために何度の引き算をしたかといえば、わずかに 1+7+3+2= 13回ということになる。これはいい。しかも2分検索と違い大きな数の自乗を計算するということがないので、16bitよりも先を計算できる。たぶん2進数でやれば31bitまでいけるのではないだろうか。

---

• 以上を2進数ベースにして簡単なプログラムにしてみた。

  bits isqrt(bits X)
  /* Xはとりあえず、 0.xxxxxxxx みたいな数値。Xにはxxxxxxxの部分が入る。 */
  /* sqrt(3)を計算したい場合は、X=.110000000000...(2進数) として、結果を2倍すればよい。
      2進数の0.11は10進数の0.75のことであり、sqrt(3/4)に相当する。これは sqrt(3)/2だから、
      結果を2倍すればいいわけだ。 */
  {
      unsigned int x = 0, y = 0, i;
      for (i = 0; i < 32; i++) {
          x = x * 4 + high(X, 2); /* xを4倍してXから上位2bitを持ってくる。
                                       次のhighで続きを持ってこられるようにXも4倍する。 */
          y *= 2;
          if (x >= y * 2 + 1) {
              x -= y * 2 + 1;
              y++;
          }
      }
      return (bits) y; /* yの前に0.を補う感じ */
  }

• これだけだと分かりにくいので、あえて2分検索版に相当する16bit版を。

  unsigned int isqrt(unsigned int xx)
  {
      unsigned int x = 0, y = 0, i;
      for (i = 0; i < 16; i++) {
          x = x << 2 | (xx >> 30); /* アセンブラならここは SHLD を使えば簡単にできる */
          xx <<= 2;
          y <<= 1;
          if (x >= y * 2 + 1) {
              x -= y * 2 + 1;
              y++;
          }
      }
      return y;
  }

• このプログラムのforのループ回数を32回にするだけで、下の桁まで求められる（答えが2^16倍になる）。

---

• これならループ回数でも2分検索版にまけてない。しかもy*2+1はy+y+1に分解できるから、これは乗算なしで開根できていることになる。いいねぇ。
• っていうか、まともに筆算アルゴリズムでやるよりもいいかも・・・というか2進数版の筆算アルゴリズムと同じになっているのかな（確認していないけど、ぱっと見てすごく似ている気はする）。

(9) 応用？

• n^3 + (3n^2+3n+1) = (n+1)^3 だから、奇数の代わりに 1, 7, 19, 37, 61,... を使うと立方根の開根にも使えそう。

こめんと欄

コメント	お名前	NameLink

   

Modified by K-tan

PukiWiki 1.4.6 Copyright © 2001-2005 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-ji
Powered by PHP 5.6.40

HTML convert time to 0.017 sec.