2010年7月10日土曜日

ビットの反転と置換

Hacker's DelightにもTAOCPにもSchroeppelのビット反転法の話が登場する. もとはHAKMEMの167項の一部にある話である.

Hacker's Delightでは(102ページ), 64ビットレジスタで7ビットxの反転を((x*0x40100401)&0x442211008)%255 と書いてあり, TAOCP(V4F1,26ページ),0<=x<2gの反転をt←((ax mod 2n)&b, y←((ct)mod 2n)>>(n-g) ただしn=g2, 0≤x<2g, a=(2n+g-1)/(2g+1-1), b=2g-1(2n-1)/(2g-1), c=(2n-g-1)/(2g-1-1) とある.

まずDelightの方.

((x*A)&B)%255
ただし A=0x40100401, B=0x442211008 と書き直す.

このAとBを二進にする.
A=100 0000 0001 0000 0000 0100 0000 0001
つまり4つの1の間に0が9個ある.
B=100 0100 0010 0010 0001 0001 0000 0000 1000

従って 7ビットgfedcbaにAを掛けると

gfedcba
1000000000100000000010000000001 A
gfedcba000gfedcba000gfedcba000gfedcba A*x
10001000010001000010001000000001000 Bでマスクすると
e a f b g c d これらが残る
43210 6543210 6543210 6543210 6543210 8ビット内の位置
11111111 28-1で割った剰余
abcdefg 反転出来た.

ビットに右から番号をつけるとaの位置は0, bは1, ..., gは6.
Aを掛けるとaの位置は0,10,20,30になり, 8の剰余は0,2,4,6になる. Bで6のaを取る. bの位置の8の剰余は1,3,5,7で5のものを取る. そういう仕掛けである.

TAOCPの方は, 定数に(2p-1)/(2q-1)が多いが, これはqビットごとの1のあるパターンを作る常套手段である.


g=3とすると, n=9. 従って

a=100010001
b=100100100
c=10101

x=rqp
これにa=100010001を掛けると
rqp0rqp0rqp
2nでmodをとると
p0rqp0rqp
b=100100100
でマスクする
p00q00r00
c=10101を掛けると
p00q00r00
p00q00r00
p00q00r00
p0pqpqrqr0r00
2nでmodをとると
pqrqr0r00
n-gビット右シフトすると
pqr
と反転出来る.


これは, 反転すべきパターンを位相をずらして複製し, 必要な部分を取り出し, Delightのように2n-1で割った剰余で揃えるか, TAOCPのように(割り算は出来ない)cをかけて混ぜ合わせ, 途中に出来た部分を取り出すかである.


TAOCPのやり方で, もっと一般的な置換が出来ないか考えたのが, 今回のブログのテーマである. 例として64ビットレジスタで8ビットのパターンを任意に置換する.

x=76543210 (ビットの名前)
とする.
置換は (0,1,2,3,4,5,6,7)->(3,2,4,1,6,0,5,7)
つまり75061423にしたい.


次のようにする.
a=0x8040201008040201
b=0xbfdfeff7fbfdfeff
y=(xa % 2^64) & b
c=0x0101010101010101
d=0x4020100804020100
z= ((xc % 2^64) >> 1) & d
m=0x14012000000a4080
(((y|z & m) * c) & (2^64-1)) >> 56

Schemeでは
(define (genperm x m) ;m permutation mask
(let* ((a #x8040201008040201) (b #xbfdfeff7fbfdfeff)
(c #x0101010101010101) (d #x4020100804020100)
(y (band (modulo (* x a) (expt 2 64)) b))
(z (band (>> (* x c) 1) d)))
(>> (band (* (band (bor y z) m) c) (- (expt 2 64) 1)) 56)))
(define m #x14012000000a4080)
実行してみると
(genperm #b11110000 m) => #b11010100
(genperm #b11001100 m) => #b10010011
(genperm #b10101010 m) => #b11001001

yはこうなる.

y= mask
76543210 ff 64ビットの内 最右の8ビット
6543210x fe
543210x7 fd
43210x76 fb
3210x765 f7
210x7654 ef
10x76543 df
0x765432 bf 最左の8ビット x: ドントケアビット

xc>>1 はこうなる.

07654321
07654321
07654321
07654321
07654321
07654321
07654321
?7654321 ?: 右シフトで左から挿入されたビット

z= ((xc mod 2^64) >>1) & d

z= mask
xxxxxxxx 00
xxxxxxx1 01
xxxxxx2x 02
xxxxx3xx 04
xxxx4xxx 08
xxx5xxxx 10
xx6xxxxx 20
x7xxxxxx 40

y|z は各列に0から7を1つずつ含む.

y|z & 0x14012000000a4080 <= 置換用マスク

76543210 7 80
65432101 5 40
54321027 1 2 0a
43210376 00
32104765 00
21057654 0 20
10676543 3 01
07765432 6 4 14
75061423

0x0101010101010101を掛ける.

7
75
75 1 2
75 1 2
75 1 2
750 1 2
75061 23
75061423

modulo 2^64 >> 56 => 75061423

完成!

置換用マスクの作り方.

00 09 18 27 36 45 54 63
08 01 10 19 28 37 46 55
56 17 02 11 20 29 38 47
48 57 26 03 12 21 30 39
40 49 58 35 04 13 22 31
32 41 50 59 44 05 14 23
24 33 42 51 60 53 06 15
16 25 34 43 52 61 62 07

y|zの行列を眺めながら作った上の表は, 0を0へ移動するには0にマスクを置く. 0を1へ移動するには9にマスクを置く. ... 0を7へ移動するには63にマスクを置く. ...7を7へ移動するには7にマスクを置く. のように読む

この表は pをqに移動するマスク位置を計算するプログラムでも計算出来る.

(define (pq p q)
(cond ((< p (+ q 1)) (+ (* -8 p) (* 9 q)))
((= p (+ q 1)) (+ (* 8 p) q))
((> p (+ q 1)) (+ 72 (* -8 p) (* 9 q)))))

上の例 0->5 1->3 2->1 3->0 4->2 5->6 6->4 7->7

(number->string
(apply +
(map (lambda (p q) (expt 2 (pq p q)))
;from pos to pos
'(0 1 2 3 4 5 6 7) '(5 3 1 0 2 6 4 7))
) 16)
=> "14012000000A4080"

0 件のコメント: