平方根なんてビット演算ですればいいじゃない

まぁ何も言わずに次のコードを実行してみましょう...

#include <stdio.h>
#include <math.h>
int mysqrt(unsigned int);

int main(){
  int i, arr[] = {0, 15, 53, 193, 456, 46340};
  for(i = 0; i < 6; ++i){
    unsigned int temp = pow(arr[i], 2);
    printf("%d = sqrt(%d) = %d\n", arr[i], temp, mysqrt(temp));
  }
  return 0;
}

int mysqrt(unsigned int x){
  int a = 0, c = 0, y = 0, i = 0, t = x;
  while(t >>= 1){
    ++i;
  }
  for(i += i & 1; i >= 0; i -= 2){
    c = (y << 1 | 1) <= x >> i;
    a = a << 1 | c;
    y = y << 1 | c;
    x -= c * y << i;
    /* if(c){
      x -= y << i;
    } */
    y += c;
  }
  return a;
}

結果はこんな感じ.

0 = sqrt(0) = 0
15 = sqrt(225) = 15
53 = sqrt(2809) = 53
193 = sqrt(37249) = 193
456 = sqrt(207936) = 456
46340 = sqrt(2147395600) = 46340

ちょっと待て...平方根がちゃんと計算出来てるじゃないか!!!
何このコード...

これ, 開平法の二進数バージョンです.
十進数の開平法については他でもいっぱい説明されているのでここでは説明しません.
開平法を行っている例(32982^2=1087812324)はこんな感じです.
binsqrt10.png
...でもこれって基数に依らずに計算できるよね ^p^
というわけで, (11011₍₂₎^2=27₍₁₀₎^2=729₍₁₀₎=1011011001₍₂₎)を二進数で開いているのがこんな感じ
binsqrt2.png
一応文章にしてみますが, 手を動かしてみないとよく分からないでしょう...

  • 元の数をxとする. これを(二進数で)二桁ずつに区切る.
  • xの上の二桁(一桁かもしれないが)を見る. 二進数なら必ず1が立つ.
  • 左に縦に1を二つ書き, その下にそれらの和10を書く.
  • 1 * 1 = 1をxの下に書き, 引き算する. 次の二桁(例では11)を下ろしてくる.
  • 10● * ● < 111となるような最大の●を見つける, というか, 二進数なら0か1なので, 大小比較.
  • 101+1=110を左の下に書き, 101 * 1を右の所に書く. 左は足し算, 右は引き算(111-101=10).
  • 次の桁(01)を下ろしてくる.
  • 今度は, 1101 > 1001なので, 0がたつ.
  • 繰り返し.
二進数で計算する利点は大小比較で立つ数字が判断できるということでしょうね.
ここのところが十進数だと結構難しいかも.
このプログラムはなんのチェックもしてないので, 負の数が入ってきたり, 46341^2 (> signed int max) が入ってきたら止まりません.

パフォーマンス? どうやって計測したらいいんだろう?
math.hのsqrtはdoubleだからこれで比較してもしょうがないし.