CRC で扱う数
CRC では、馴染みのない性質を持つ数を使って計算します。
2値循環数
整数の並びは
\cdots,-2,-1,0,+1,+2,\cdots
ですが、2 で割ったときの余りに置き換えると
\cdots,\{-2→0\},\{-1→1\},\{0→0\},\{+1→1\},\{+2→0\},\cdots
0 と 1 の並び
\cdots,0,1,0,1,0,\cdots
になります。これを「 2値循環数 」と呼ぶことにします。整数 $N$ に対する、2値循環数 $n$ は
n \equiv N \mod 2
です。2値循環数では、偶数(0)または奇数(1)以外を考慮しません。
2値循環数の加法
2つの 2値循環数 $a,b$ を
\displaylines {
a \equiv A \mod 2 \\
b \equiv B \mod 2
}
とすると、2値循環数の世界での $a$ と $b$ の足し算{ $c=a+b$ }は
c = a + b \equiv \left( A + B \right) \equiv (a + b) \mod 2
になります。引き算{ $d=a-b$ }は
d = a - b \equiv (A - B) \equiv (a - b) \mod 2
になります。$a$ と $b$ はそれぞれ 2 通りだから、計算は4通り
| a | b | c | d |
|---|---|---|---|
| 0 | 0 | 0 | 0 |
| 0 | 1 | 1 | 1 |
| 1 | 0 | 1 | 1 |
| 1 | 1 | 0 | 0 |
で、足し算(c)と引き算(d)は、論理演算の排他的論理和(XOR)と同じ計算になります。
整数の計算と違うところは
\displaylines {
1 + 1 \equiv 0 \equiv +2 \mod 2 \\
0 - 1 \equiv 1 \equiv -1 \mod 2
}
ですが、{0=偶数,1=奇数} で考えると
| $A$ | $B$ | $A+B$ | $A-B$ |
|---|---|---|---|
| 偶数 | 偶数 | 偶数 | 偶数 |
| 偶数 | 奇数 | 奇数 | 奇数 |
| 奇数 | 偶数 | 奇数 | 奇数 |
| 奇数 | 奇数 | 偶数 | 偶数 |
になります。整数で計算しても、それが「偶数か?奇数か?」に注目します。
2値循環数を係数とする多項式
変数 $x$ と係数 $a_n$ の多項式を
A = a_0 x^0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots + a_k x^k = \sum_{n=0}^k a_n x^n
とします。ここで、係数 $a_n$ は2値循環数とします。すると、係数は { 0 または 1 } に限られます。さらに、係数 {$a_0,a_1,\cdots,a_k$} を逆順に並べて
A_{(2)} = a_k \cdots a_3 a_2 a_1 a_0
2進数表現にします。( $x=2$ を $A$ に代入すると $A=A_{(2)}$ になります)
「2値循環数を係数とする多項式」の計算を、2進数表現に簡略化して行います。
足し算
2値循環数を係数にした多項式 $A,B$ を
\displaylines {
A = a_0 x^0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots + a_k x^k = \sum_{n=0}^k a_n x^n \\
B = b_0 x^0 + b_1 x^1 + b_2 x^2 + b_3 x^3 + \cdots + b_k x^k = \sum_{n=0}^k b_n x^n
}
として、$A$ と $B$ の足した { $C=A+B$ } は
\displaylines {
C = A + B = \sum_{n=0}^k \left( a_n + b_n \right) x^n = \sum_{n=0}^k c_n x^n \\
\{ c_n = a_n + b_n = XOR \left( a_n, b_n \right) \}
}
になります。係数を並べた2進数表現では
\displaylines {
A_{(2)} = a_k \cdots a_3 a_2 a_1 a_0 \\
B_{(2)} = b_k \cdots b_3 b_2 b_1 b_0 \\
C_{(2)} = c_k \cdots c_3 c_2 c_1 c_0 \\
C_{(2)} = XOR \left( A_{(2)}, B_{(2)} \right)
}
となるので、排他的論理和(XOR)で計算できます。
引き算
足し算と同様に、引き算 { $C=A-B$ } は
\displaylines {
C = A - B = \sum_{n=0}^k \left( a_n - b_n \right) x^n = \sum_{n=0}^k c_n x^n \\
\{ c_n = a_n - b_n = XOR \left( a_n, b_n \right) \}
}
となって、 足し算と同じ計算
C_{(2)} = XOR \left( A_{(2)}, B_{(2)} \right)
になります。
掛け算
$A$ に $B$ を掛けた { $C=A\ B$ } は
\displaylines {
C = A\ B = \left( \sum_{n=0}^k a_n x^n \right) \left( \sum_{n=0}^k b_n x^n \right) \\
= \left( a_0 x^0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots + a_k x^k \right) b_0 x^0 \\
+ \left( a_0 x^0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots + a_k x^k \right) b_1 x^1 \\
+ \left( a_0 x^0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots + a_k x^k \right) b_2 x^2 \\
+ \left( a_0 x^0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots + a_k x^k \right) b_3 x^3 \\
+ \cdots \\
+ \left( a_0 x^0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots + a_k x^k \right) b_k x^k
}
となりますが、$x$ の指数は、2進数表現の桁位置に相当することを念頭に入れつつ
\displaylines {
= \left( a_0 x^0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots + a_n x^{k+0} \right) b_0 \\
+ \left( a_0 x^1 + a_1 x^2 + a_2 x^3 + a_3 x^4 + \cdots + a_n x^{k+1} \right) b_1 \\
+ \left( a_0 x^2 + a_1 x^3 + a_2 x^4 + a_3 x^5 + \cdots + a_n x^{k+2} \right) b_2 \\
+ \left( a_0 x^3 + a_1 x^4 + a_2 x^5 + a_3 x^6 + \cdots + a_n x^{k+3} \right) b_3 \\
+ \cdots \\
+ \left( a_0 x^k + a_1 x^{k+1} + a_2 x^{k+2} + a_3 x^{k+3} + \cdots + a_n x^{2k} \right) b_k
}
とします。2進数表現にすると
\displaylines {
C_{(2)} = A_{(2)} \times B_{(2)} \\
= \left( a_k \cdots a_3a_2a_1a_0 \right) \times b_0 \\
+ \left( a_k \cdots a_3a_2a_1a_0\ 0 \right) \times b_1 \\
+ \left( a_k \cdots a_3a_2a_1a_0\ 00 \right) \times b_2 \\
+ \left( a_k \cdots a_3a_2a_1a_0\ 000 \right) \times b_3 \\
+ \cdots \\
+ \left( a_k \cdots a_3a_2a_1a_0\ 0_k \cdots 0_1 \right) \times b_k
}
だから、整数の筆算と同じように計算できます。しかし、ここでの足し算は XOR です。
例: 1101 * 0101 = 0111001
1101 : 1101 * 0001 1101 : 1101 * 1
0000 : 1101 * 0000 00000 : 11010 * 0
1101 : 1101 * 0100 110100 : 110100 * 1
0000 : 1101 * 0000 0000000 : 1101000 * 0
--------:----------- --------:------------
0111001 : 1101 * 0101 0111001 : 1101 * 0101
割り算
$A$ を $B$ で割ったときの { 商, 余 } を { $C,D$ } とすると
\displaylines {
D \equiv A \mod B \\
C = \frac{A-D}{B}
}
ですが、多項式の計算では、$A$ から $D$ を引いたものを因数分解するイメージ
A - D = B\ C
になります。しかし、2進数表現における掛け算の計算手順を考えると、割り算も整数の筆算と同じように計算できることになります。ただし、引き算は XOR です。
例: 111010 ÷ 101 = 1101 余り 11
1101
--------
101 ) 111010
101
----
100
101
----
110
101
---
11
整数の計算と違って、「除数と同じ桁数の数」になった時点で、引き算(XOR)ができます。
割り算の余り
割り算の余りは、除数の桁数と同じ桁数にはなりません。また、余りの範囲は、除数の桁数で決まり、 同じ桁数で、異なる除数でも、余りの範囲は同じ という特徴があります。よって、 余りの範囲を決めてから、除数を決める ことができます。
除数 : 余り
1x : 1
1xx : 01, 10, 11
1xxx : 001, 010, 011, 100, 101, 110, 111
...
除数の選択
因数分解できる多項式は、2進数表現では掛け算になり、整数の合成数に相当します。となると、整数の素数に相当する多項式がありそうです。それは、既約多項式とよばれています。
様々な数を割って、余りの分布に偏りが少ないものを除数として選択します。整数では素数になるので、既約多項式を選択することになります。簡単な既約多項式の求め方として、$A \times B$ で合成数(相当)を排除した残り(エラトステネスのふるい相当)があります。
CRC の計算
除数を $B$ とし、$B$ の桁数を $k+1$ とすると、余りの範囲は $k$ 桁になります。データ列 $D$ ($m$桁)に $k$桁のゼロを追加して、被除数 $A$ とします。そして剰余 { $C \equiv A \mod B$ } を求めます。
\displaylines {
B_{(2)} = b_k b_{k-1} \cdots b_2 b_1 b_0 \\
D_{(2)} = d_{m-1} \cdots d_2 d_1 d_0 \\
A_{(2)} = d_{m-1} \cdots d_2 d_1 d_0 \ 0_k \cdots 0_1 \\
C \equiv A \mod B \\
C_{(2)} = c_{k-1} c_{k-2} \cdots c_2 c_1 c_0
}
この計算が基本形です。
B = bb..bb
D = dd....dd
A = dd....dd 00..00
C = cc..cc = A % B
データ(バイト)列の並び順
バイト配列 data[] をビット位置 bit_index で素直にアクセスすると
bit_access_el.c
int bit_access_el(const unsigned char *data, size_t bit_index)
{
size_t pos_byte = bit_index >> 3;
size_t pos_bit = bit_index & 7;
return (data[pos_byte] >> pos_bit) & 1;
}
になります。よって、データ(バイト)列は下位の桁から並んでいるように見えます。
データ例
1 バイト目のデータ : 10100001 = 0xa1
ビット 0 = 1
ビット 1 = 0
ビット 2 = 0
ビット 3 = 0
ビット 4 = 0
ビット 5 = 1
ビット 6 = 0
ビット 7 = 1
2 バイト目のデータ : 10011010 = 0x9a
ビット 0 = 0
ビット 1 = 1
ビット 2 = 0
ビット 3 = 1
ビット 4 = 1
ビット 5 = 0
ビット 6 = 0
ビット 7 = 1
...
バイト並び 0xa1,0x9a,... が、大きな整数 0x...9aa1 に見える(リトルエンディアン)
割り算は上位の桁から計算しますが、データは下位の桁から並んでいて、都合がよくありません。そこで、ビット並びを修正して
bit_access_eb.c
int bit_access_eb(const unsigned char *data, size_t bit_index)
{
size_t pos_byte = bit_index >> 3;
size_t pos_bit = 7 - (bit_index & 7); /* ビット並びを修正 */
return (data[pos_byte] >> pos_bit) & 1;
}
とすると
バイト並び 0xa1,0x9a,... が、大きな整数 0xa19a... に見える(ビッグエンディアン)
データは上位の桁から並ぶように見えます。このバイト列をデータ列として計算する方法を左算とします。
「桁の並び順は逆」とする
「左の桁を上位」とする概念を横に置いて、「 右の桁を上位 」とします。これを右算とします。
「足し算(XOR)」と「引き算(XOR)」は計算が(多項式の係数だから)桁ごとに完結しているので、桁の並び順は関係ありません。「掛け算」と「割り算」は、筆算の方法を左右反転します。
例を「左桁上位」と「右桁上位」の併記にします。
掛け算
例: 1101 * 101 = 111001
ビット列の反転操作
bit_reflect( 1101) = 1011
bit_reflect( 101) = 101
bit_reflect(111001) = 100111
左桁上位 右桁上位
1101 : 101 * 001 100 * 101 : 1011
0000 : 101 * 000 000 * 101 : 0000
1101 : 101 * 100 001 * 101 : 1011
-------:----------- -----------:-------
111001 : 1101 * 101 101 * 1011 : 100111
割り算
例: 111011 ÷ 110 = 1011 余り 1
ビット列の反転操作
bit_reflect(111011) = 110111
bit_reflect( 110) = 011
bit_reflect( 1011) = 1101
bit_reflect( 1) = 1
左桁上位 右桁上位
1011 1101
-------- -------
110 ) 111011 110111 ( 011
110 011
----- -----
101 101
110 011
---- ----
111 111
110 011
--- ---
1 1
右算では
バイト並び 0xa1,0x9a,... が、大きな整数 bit_reflect(0x...9aa1) に見える
となるので、データは上位の桁から並ぶように見えます。
左算と右算では、結果が違います。
CRC 計算(基本形)のプログラムを書いてみる
開発環境: Linux [Debian 9 (x64)]
言語: C
コンパイラ: GCC 6.3.0
除数が $2^{32}$ 未満のプログラム
crc.c
#include <stdio.h>
#include <string.h>
typedef unsigned int uint32;
/* ビット列の反転
* d: データ
* b: ビット数
*/
uint32 bit_reflect(uint32 d, int b)
{
uint32 r = 0;
while (b-- > 0)
{
r <<= 1;
r |= d & 1;
d >>= 1;
}
return r;
}
/* 最上位ビットの位置を得る
* x: 調べる値
*/
int bsr(uint32 x)
{
#if 0
int b;
for (b = 31; b >= 0; b--)
if (((x >> b) & 1) != 0) break;
return b;
#else /* IA-32 では上の検索処理を bsr 命令で行える */
int b = -1;
__asm__ ("bsr %1,%0": "+r" (b) : "r" (x));
return b;
#endif
}
/* 除算(1ビット分)
* r: 計算中の剰余
* p: 除数
* d: データ(1ビット)
*/
uint32 crc_step1(uint32 r, uint32 p, uint32 d)
{
int k = bsr(p); /* 余りの桁数 */
r = (r << 1) | (d & 1); /* 計算中の剰余を桁送り */
if ((r & (1 << k)) != 0) /* 引き算(XOR)が可能か？ */
r ^= p; /* 引き算(XOR) */
return r; /* 更新された剰余 */
}
/* 除算(8ビット分：左算用)
* r: 計算中の剰余
* p: 除数
* d: データ(8ビット)
*/
uint32 crc_step8_left(uint32 r, uint32 p, uint32 d)
{
int b;
for (b = 0; b < 8; b++, d <<= 1)
r = crc_step1(r, p, (d >> 7)); /* d: bit 7 を bit 0 へ移動 */
return r;
}
/* 除算(8ビット分：右算用)
* r: 計算中の剰余
* p: 除数
* d: データ(8ビット)
*/
uint32 crc_step8_right(uint32 r, uint32 p, uint32 d)
{
int b;
for (b = 0; b < 8; b++, d >>= 1)
r = crc_step1(r, p, d);
return r;
}
/* CRC 左算
* s: データ列
* n: バイト数
* p: 除数
*/
uint32 crc_left(const char *s, size_t n, uint32 p)
{
int k = bsr(p), b = k; /* 余りの桁数 */
uint32 r = 0; /* 剰余: 初期値は 0 */
while (n-- > 0) r = crc_step8_left(r, p, *s++); /* データ列 D */
while (b-- > 0) r = crc_step1(r, p, 0); /* 0 の列を追加して A 化 */
return r; /* 剰余 */
}
/* CRC 右算
* s: データ列
* n: バイト数
* p: 除数
*/
uint32 crc_right(const char *s, size_t n, uint32 p)
{
int k = bsr(p), b = k; /* 余りの桁数 */
uint32 r = 0; /* 剰余: 初期値は 0 */
while (n-- > 0) r = crc_step8_right(r, p, *s++); /* データ列 D */
while (b-- > 0) r = crc_step1(r, p, 0); /* 0 の列を追加して A 化 */
return r; /* 剰余 */
}
/* CRC テスト
* s: 文字列
* p: 除数(2^32 未満)
*/
void crc_test(const char *name, const char *s, uint32 p)
{
int k = bsr(p); /* 余りの桁数 */
size_t n = strlen(s); /* 文字列の長さ */
printf("**** %s: 文字列[%zd]='%s' \n", name, n, s);
printf("CRC 左算: 除数=%#010x, 剰余=%#010x\n", p, crc_left(s, n, p));
printf("CRC 右算: 除数=%#010x, 剰余=%#010x\n", bit_reflect(p, k), bit_reflect(crc_right(s, n, p), k));
}
/* メイン関数 */
int main(int argc, char **argv)
{
const char *data;
if (argc != 2)
{
printf("Usage: %s data\n", argv[0]);
return 1;
}
data = argv[1];
crc_test("CRC-8", data, 0x00131);
crc_test("CRC-16", data, 0x11021);
crc_test("CRC-30", data, 0x30185CE3);
return 0;
}
ビルドと実行結果
$ cc -W -Wall -O2 -o crc crc.c
$ ./crc 'test:abcde'
**** CRC-8: 文字列[10]='test:abcde'
CRC 左算: 除数=0x00000131, 剰余=0x000000f9
CRC 右算: 除数=0x0000008c, 剰余=0x0000001b
**** CRC-16: 文字列[10]='test:abcde'
CRC 左算: 除数=0x00011021, 剰余=0x000057ec
CRC 右算: 除数=0x00008408, 剰余=0x0000b8dd
**** CRC-30: 文字列[10]='test:abcde'
CRC 左算: 除数=0x30185ce3, 剰余=0x1549a57f
CRC 右算: 除数=0x18e74301, 剰余=0x17fca10b
除数の範囲が $[ 2^{32}, 2^{33} )$ のプログラム (CRC32命令を使用)
crc32.c
#include <stdio.h>
#include <string.h>
typedef unsigned int uint32;
/* ビット列の反転
* d: データ
* b: ビット数
*/
uint32 bit_reflect(uint32 d, int b)
{
uint32 r = 0;
while (b-- > 0)
{
r <<= 1;
r |= d & 1;
d >>= 1;
}
return r;
}
/* 除算(1ビット分)
* r: 計算中の剰余
* p: 除数の下位32ビット
* d: データ(1ビット)
*/
uint32 crc32_step1(uint32 r, uint32 p, uint32 d)
{
uint32 u = (r >> 31);
r = (r << 1) | (d & 1);
if (u != 0) r ^= p;
return r; /* 更新された剰余 */
}
/* 除算(8ビット分：左算用)
* r: 計算中の剰余
* p: 除数の下位32ビット
* d: データ(8ビット)
*/
uint32 crc32_step8_left(uint32 r, uint32 p, uint32 d)
{
int b;
for (b = 0; b < 8; b++, d <<= 1)
r = crc32_step1(r, p, (d >> 7)); /* d: bit 7 を bit 0 へ移動 */
return r;
}
/* 除算(8ビット分：右算用)
* r: 計算中の剰余
* p: 除数の下位32ビット
* d: データ(8ビット)
*/
uint32 crc32_step8_right(uint32 r, uint32 p, uint32 d)
{
int b;
for (b = 0; b < 8; b++, d >>= 1)
r = crc32_step1(r, p, d);
return r;
}
/* CRC 左算
* s: データ列
* n: バイト数
* p: 除数の下位32ビット
*/
uint32 crc32_left(const char *s, size_t n, uint32 p)
{
int b = 32;
uint32 r = 0; /* 剰余: 初期値は 0 */
while (n-- > 0) r = crc32_step8_left(r, p, *s++); /* データ列 D */
while (b-- > 0) r = crc32_step1(r, p, 0); /* 0 の列を追加して A 化 */
return r; /* 剰余 */
}
/* CRC 右算
* s: データ列
* n: バイト数
* p: 除数
*/
uint32 crc32_right(const char *s, size_t n, uint32 p)
{
int b = 32;
uint32 r = 0; /* 剰余: 初期値は 0 */
while (n-- > 0) r = crc32_step8_right(r, p, *s++); /* データ列 D */
while (b-- > 0) r = crc32_step1(r, p, 0); /* 0 の列を追加して A 化 */
return r; /* 剰余 */
}
/* CRC32 命令による右算
* s: データ列
* n: バイト数
* 除数は 0x11edc6f41
*/
uint32 x86_crc32(const char *s, size_t n)
{
uint32 r = 0;
while (n-- > 0)
__asm__ ("crc32b %1,%0": "+r" (r) : "g" (*s++));
return r;
}
/* CRC テスト
* s: 文字列
* p: 除数 [2^32,2^33) の下位32ビット
*/
void crc32_test(const char *name, const char *s, uint32 p)
{
size_t n = strlen(s); /* 文字列の長さ */
printf("**** %s: 文字列[%zd]='%s' \n", name, n, s);
printf("左算: 除数=0x1%08x, 剰余=%#010x\n", p, crc32_left(s, n, p));
printf("右算: 除数=0x1%08x, 剰余=%#010x\n", bit_reflect(p, 32), bit_reflect(crc32_right(s, n, p), 32));
}
/* CRC-32C テスト
* s: 文字列
*/
void x86_crc32_test(const char *s)
{
size_t n = strlen(s); /* 文字列の長さ */
printf("**** CRC-32C 右算: 文字列[%zd]='%s' \n", n, s);
printf("通常計算: 剰余=%#010x\n", bit_reflect(crc32_right(s, n, 0x1edc6f41), 32));
printf("専用計算: 剰余=%#010x\n", x86_crc32(s, n));
}
/* SSE4.2 の調査(簡易版) */
int check_sse42()
{
uint32 ecx;
__asm__
("mov $1,%%eax;"
"cpuid;"
"mov %%ecx,%0"
: "=r" (ecx)
: /**/
: "eax", "ebx", "ecx", "edx");
return ((ecx >> 20) & 1);
}
/* メイン関数 */
int main(int argc, char **argv)
{
const char *data;
if (argc != 2)
{
printf("Usage: %s data\n", argv[0]);
return 1;
}
data = argv[1];
crc32_test("CRC-32", data, 0x04c11db7);
crc32_test("CRC-32K", data, 0x741b8cd7);
if (check_sse42()) x86_crc32_test(data);
else printf("SSE4.2: CRC32 命令は使用できません\n");
}
ビルドと実行結果
$ cc -W -Wall -O2 -o crc32 crc32.c
$ ./crc32 'test:abcde'
**** CRC-32: 文字列[10]='test:abcde'
左算: 除数=0x104c11db7, 剰余=0xa95c1a25
右算: 除数=0x1edb88320, 剰余=0x35ea88f1
**** CRC-32K: 文字列[10]='test:abcde'
左算: 除数=0x1741b8cd7, 剰余=0x1bdbe7eb
右算: 除数=0x1eb31d82e, 剰余=0xa73e4a2c
**** CRC-32C 右算: 文字列[10]='test:abcde'
通常計算: 剰余=0xd5416383
専用計算: 剰余=0xd5416383
となったので、組んだプログラムと CRC32 命令の計算は合っているようです。
SSE4.2 拡張の CRC32 命令は、除数が固定 ( 0x11EDC6F41 ) だから、除数が違う CRC32 ( 0x104C11DB7) { gzip, zlib, PNG など } には適用できないようです。
CRC の値で CRC 計算する (M系列乱数?)
CRC は剰余なので、フェルマーの小定理
a^{p-1} \equiv 1 \mod p
からの連想で、数列の循環が気になります。
初期値を $x_0 = 1$ として、再帰的に CRC32 命令を実行
x_{n+1} = CRC32C \left( x_n \right)
したら、$x_0 = x_\left( 2^{31} - 1 \right)$ になりました。このとき、使った関数は
uint32 CRC32C(uint32 x)
{
uint32 y = 0;
__asm__ ("crc32l %1,%0" : "+r" (y) : "r" (x));
return y;
}
です。ここで、$y$ を弄って
uint32 CRC32C(uint32 x)
{
uint32 y = 1; /* 変更 */
__asm__ ("crc32l %1,%0" : "+r" (y) : "r" (x));
return y;
}
とすると、$x_0 = x_\left( 2^{32} - 2 \right)$ となって、uint32 で出現しない値は
0x7181278e
0x8d250a21
の 2つでした。$y$ の値によって、出現しない値が変わります。
これなら、簡単に作れる乱数として使えそうです。