今朝書いたNLZのコードについてのメモ。

朝見たときは全然意味不明コードで内容を理解しないままNLZ*1に適用しましたが今日1日でそれなりに分かりました。

１．"\0\1\2\xf..."って何？

多分0〜31の配列に相当するんだろうなぁとは思っていましたが \xf なんて表記、初めて見たので何よそれ？って感じでしたがやっぱり予想通りでした。

なので

unsigned nlzM(unsigned x) {
   if (x == 0) { return 32; }
   x = x | (x >> 1);
   x = x | (x >> 2);
   x = x | (x >> 4);
   x = x | (x >> 8);
   x = x | (x >>16);
   return 31 - ("\0\1\2\xf\x1d\3\x17\x10\x1e\x1b\4\6\xc\x18\x8\x11"
                "\x1f\xe\x1c\x16\x1a\5\xb\7\xd\x15\x19\xa\x14\x9\x13\x12")
               [0x5763e69U * (x - (x >> 1)) >> 27];
}

を

unsigned nlzM(unsigned x) {
   static const unsigned char tbl[32] =
   { 0,1,2,15,29,3,23,16,30,27,4,6,12,24,8,17,31,14,28,22,26,5,11,7,13,21,25,10,20,9,19,18 };

   if (x == 0) { return 32; }
   x = x | (x >> 1);
   x = x | (x >> 2);
   x = x | (x >> 4);
   x = x | (x >> 8);
   x = x | (x >>16);
    return 31 - tbl[0x5763e69U * x >> 27];
}

と書き換えても得られる結果は同じです。

ただし

unsigned nlzM(unsigned x) {
   static const unsigned tbl[32] =
   { 0,1,2,15,29,3,23,16,30,27,4,6,12,24,8,17,31,14,28,22,26,5,11,7,13,21,25,10,20,9,19,18 };

   if (x == 0) { return 32; }
   x = x | (x >> 1);
   x = x | (x >> 2);
   x = x | (x >> 4);
   x = x | (x >> 8);
   x = x | (x >>16);
    return 31 - tbl[0x5763e69U * x >> 27];
}

とすれば若干速くなります。（サイズは増えるけど）

２．0x5763e69はどうやって求めた？？

M系列というキーワードがあちらこちらで見られるので、まあコイツを使えば出て来るんでしょうが、その実際はってことで色々試した結果のコードがコチラ！

void MagicNum(void)
{
    unsigned i, j, flg, val;
    unsigned arr[32];
    char str[4*32+1];
    
    val=0;
    arr[0]=0;
    flg=1;
    
    for (i=1; i<32; i++)
    {
        val = (val<<1) | flg;
        arr[i] = val & 0x1f;
        flg = 1 & ((val>>4) + (val>>1));                       /* no.1 (5,2) */
        //flg = 1 & ((val>>4) + (val>>3) + (val>>2) + (val>>1)); /* no.2 (5,4,3,2) */
        //flg = 1 & ((val>>4) + (val>>3) + (val>>1) + (val>>0)); /* no.3 (5,4,2,1) */
        //flg = 1 & ((val>>4) + (val>>2));                       /* no.4 (5,3) */
    }
    val>>=4;
    sprintf(str, "\\0");
    for (i=1; i<32; i++)
    {
        for (j=1; arr[j]!=i; j++){}
        if (j>7) { sprintf(str, "%s\\x%x", str, j); }
        else     { sprintf(str, "%s\\%x" , str, j); }
    }
    
    printf("hex : 0x%x\n", val);
    printf("str : %s"    , str);
    getchar();
}

ここで使われるのは5次のM系列らしいです。
ではその5次のM系列ってどうすんのよ？ってのはM-series generator polynomialsを参考にしました。

上記コードの(5,2)ってところの意味は、今回値を1として5→4個前の0/1の値と2→1個前の0/1の値のXorが次回の0/1ってことです。
この(5,2)の式で0x5763e69が求められます。

なので同様に(5,4,3,2)なら

unsigned nlzM(unsigned x) {
   if (x == 0) { return 32; }
   x = x | (x >> 1);
   x = x | (x >> 2);
   x = x | (x >> 4);
   x = x | (x >> 8);
   x = x | (x >>16);
   return 31 - ("\0\1\2\xd\x19\3\x1c\xe\xb\x1a\x9\4\x1d\6\xf\x13\x1f\xc\x18\x1b\xa\x8\5\x12\x1e\x17\7\x11\x16\x10\x15\x14")
               [0x5a8ef93U * l_tmp_ul >> 27];
}

(5,4,2,1)では

   return 31 - ("\0\1\xe\2\xb\xf\3\x1b\xc\x9\7\x10\x18\4\x1c\x12\x1f\xd\xa\x1a\x8\6\x17\x11\x1e\x19\5\x16\x1d\x15\x14\x13")
               [0x6a45f67U * l_tmp_ul >> 27];

というコードが得られます。

ちなみにM-series generator polynomialsでは上記の3パターン以外でも、(5,3)で

   return 31 - ("\0\1\2\x13\3\6\x14\xc\x1e\4\x1c\7\x9\x15\x18\xd\x1f\x12\5\xb\x1d\x1b\x8\x17\x11\xa\x1a\x16\x10\x19\xf\xe")
               [0x4b3e375U * l_tmp_ul >> 27];

NLZ的にはOKみたいです。（M系列的に正しい組み合わせかは不明）
ちなみに上記4パターンでの最終的なNLZの結果・処理速度には何ら影響はしないので、まあどれでもいいと思われます。

改めてM系列をこれに適用しようという発想が素晴らしいなと実感しました。

*1:Number of Leading Zero