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

バックアップ一覧
差分を表示
現在との差分を表示
ソースを表示
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 isqrt(int x)
{
    int a = 0, b = 46340, c;
    while (a + 1 < b) {
        c = (a + b) / 2;
        if (x < c * c)
            b = c;
        else
            a = c;
    }
    return a;
}

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

↑

(6) 改良

（まだ書いてない）

↑

(9) 応用？

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

↑

こめんと欄

コメント	お名前	NameLink