それなりに速い自前の平方根関数(sqrt)を作りました。

仕事でたまに組み込み用（浮動小数点の使えない整数オンリーな環境）のプログラムを触るのですが、精度が要求される演算などの計算では浮動小数点を使わざるを得ないところもあり、その箇所では自前の四則演算用の関数を使って処理をしています。
平方根の算出関数はこれまでは使うことがありませんでしたが、今後のことを考え勉強がてら作ってみました。

浮動小数点について

floatの変数をunsigned longとして参照する*1と
左から符号1bit、指数部8bit、仮数部が23bit（＋隠しbitを加えて24bit）
0x3F800000（1.0f）を例にとると
符号：0
指数部：0x7F(127)
仮数部：0
となり
0x3F800000 → (-1)^0 * 2^(127-127) * (1 + 0) = 1.0になります。
（端折りまくっていますので興味のある方はIEEE 754 - Wikipediaをどうぞ）

平方根（ルート）の求め方

平方根の指数は引数の指数を2で割った値になります。
例えば √(2^10) なら → 2^5です。

ちなみに指数が奇数のときは2で割って切り捨てます。
その代わりに仮数部を2倍してから仮数部のルートを求めます。
例．√(1.9×2^11)＝√3.8×2^5

次に仮数部について
仮数部は元の値の仮数部24bitを2^8倍します。
そして指数が偶数のときは2で割ります。

・指数が偶数ならとり得る値の範囲は 0x40000000〜0x7FFFFF80
・指数が奇数ならとり得る値の範囲は 0x80000000〜0xFFFFFF00

これで前準備は完了でこの32bitの値の平方根16bitを計算します。
テンポラリの値 y に初期値2^15をセットし
32bit値 ≧ y * y なら y = y + 2^15 そうでないなら y = y - 2^15
32bit値 ≧ y * y なら y = y + 2^14 そうでないなら y = y - 2^14
：
(中略)
：
32bit値 ≧ y * y なら y = y + 2^0 そうでないなら y = y - 2^0

これで16bit分が求まりました。
浮動小数点の仮数部は24bitなので、残り8bit必要です。
この8bitは以下のようにして求めます。

32bit値 * 2^0 ≧ y * y なら y = y * 2 + 1 そうでないなら y = y * 2 - 1
32bit値 * 2^2 ≧ y * y なら y = y * 2 + 1 そうでないなら y = y * 2 - 1
：
(中略)
：
32bit値 * 2^7 ≧ y * y なら y = y * 2 + 1 そうでないなら y = y * 2 - 1

最初の16bitのところで初期値2^15に対し2^14、2^13、、、2^0を
次の8bitは1/2^1(0.5)、1/2^2(0.25)、、、1/2^7(0.0078125)を条件に従って加算／減算し
その結果に2^8を掛けているようなイメージです。
これを実際にコーディングしたものが以下のソースです。

unsigned long sqr (unsigned long x)
{
    unsigned long exponent = (x >> 23) - 127;        /*  指数部 :  8bit  */
    unsigned long mantissa = (x <<  8) | 0x80000000; /*  仮数部 : 32bit  */
    unsigned long y;
    
    
    /*  引数が0以下の場合は0を返す  */
    if ((signed)x <= 0)
    {
        return 0;
    }
    
    /*  指数部が偶数の場合、仮数を半分にする  */
    if ((exponent & 1) == 0)
    {
        mantissa >>= 1;
    }
    
    /*  指数部を半分  */
    exponent = (exponent >> 1) + 127;
    
    /*  仮数部16bit分計算  */
    y = (1 << 15) + (1 << 14);
    for (i = (1 << 13); i > 0; i >>= 1) {
        if (mantissa >= y * y) {
            y += i;
        } else {
            y -= i;
        }
    }
    
    /*  残り8+2bit分計算  */
    for (i = 0; i < 10*2; i += 2) {
        if ((signed)((mantissa << i) - y * y) >= 0) {  /*  オーバーフローするが結果は正しい  */
            y = (y << 1) + 1;
        } else {
            y = (y << 1) - 1;
        }
    }
    
    
    /*  26-2=24bit  */
    return  (exponent << 23)
          | (((y + 1) >> 2) & 0x7FFFFF);
}