-
Notifications
You must be signed in to change notification settings - Fork 0
/
bch15_5.c
184 lines (169 loc) · 5.12 KB
/
bch15_5.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
/*Copyright (C) 2008-2009 Timothy B. Terriberry (tterribe@xiph.org)
You can redistribute this library and/or modify it under the terms of the
GNU Lesser General Public License as published by the Free Software
Foundation; either version 2.1 of the License, or (at your option) any later
version.*/
#include "bch15_5.h"
/*A cycle in GF(2**4) generated by alpha=(x**4+x+1).
It is extended an extra 16 entries to avoid some expensive mod operations.*/
static const unsigned char gf16_exp[31]={
1,2,4,8,3,6,12,11,5,10,7,14,15,13,9,1,2,4,8,3,6,12,11,5,10,7,14,15,13,9,1
};
/*The location of each integer 1...16 in the cycle.*/
static const signed char gf16_log[16]={
-1,0,1,4,2,8,5,10,3,14,9,7,6,13,11,12
};
/*Multiplication in GF(2**4) using logarithms.*/
static unsigned gf16_mul(unsigned _a,unsigned _b){
return _a==0||_b==0?0:gf16_exp[gf16_log[_a]+gf16_log[_b]];
}
/*Division in GF(2**4) using logarithms.
The result when dividing by zero is undefined.*/
static unsigned gf16_div(unsigned _a,unsigned _b){
return _a==0?0:gf16_exp[gf16_log[_a]+15-gf16_log[_b]];
}
/*Multiplication in GF(2**4) when the second argument is known to be non-zero
(proven by representing it by its logarithm).*/
static unsigned gf16_hmul(unsigned _a,unsigned _logb){
return _a==0?0:gf16_exp[gf16_log[_a]+_logb];
}
/*The syndrome normally has five values, S_1 ... S_5.
We only calculate and store the odd ones in _s, since S_2=S_1**2 and
S_4=S_2**2.
Returns zero iff all the syndrome values are zero.*/
static int bch15_5_calc_syndrome(unsigned _s[3],unsigned _y){
unsigned p;
int i;
int j;
p=0;
for(i=0;i<15;i++)if(_y&1<<i)p^=gf16_exp[i];
_s[0]=p;
p=0;
for(i=0;i<3;i++)for(j=0;j<5;j++)if(_y&1<<5*i+j)p^=gf16_exp[j*3];
_s[1]=p;
p=0;
for(i=0;i<5;i++)for(j=0;j<3;j++)if(_y&1<<3*i+j)p^=gf16_exp[j*5];
_s[2]=p;
return _s[0]!=0||_s[1]!=0||_s[2]!=0;
}
/*Compute the coefficients of the error-locator polynomial.
Returns the number of errors (the degree of the polynomial).*/
static int bch15_5_calc_omega(unsigned _o[3],unsigned _s[3]){
unsigned s02;
unsigned tt;
unsigned dd;
int d;
_o[0]=_s[0];
s02=gf16_mul(_s[0],_s[0]);
dd=_s[1]^gf16_mul(_s[0],s02);
tt=_s[2]^gf16_mul(s02,_s[1]);
_o[1]=dd?gf16_div(tt,dd):0;
_o[2]=dd^gf16_mul(_s[0],_o[1]);
for(d=3;d>0&&!_o[d-1];d--);
return d;
}
/*Find the roots of the error polynomial.
Returns the number of roots found, or a negative value if the polynomial did
not have enough roots, indicating a decoding error.*/
static int bch15_5_calc_epos(unsigned _epos[3],unsigned _s[3]){
unsigned o[3];
int nerrors;
int d;
int i;
d=bch15_5_calc_omega(o,_s);
nerrors=0;
if(d==1)_epos[nerrors++]=gf16_log[o[0]];
else if(d>0){
for(i=0;i<15;i++){
int i2;
i2=gf16_log[gf16_exp[i<<1]];
if(!(gf16_exp[i+i2]^gf16_hmul(o[0],i2)^gf16_hmul(o[1],i)^o[2])){
_epos[nerrors++]=i;
}
}
if(nerrors<d)return -1;
}
return nerrors;
}
int bch15_5_correct(unsigned *_y){
unsigned s[3];
unsigned epos[3];
unsigned y;
int nerrors;
int i;
y=*_y;
if(!bch15_5_calc_syndrome(s,y))return 0;
nerrors=bch15_5_calc_epos(epos,s);
if(nerrors>0){
/*If we had a non-zero syndrome value, we should always find at least one
error location, or we've got a decoding error.*/
for(i=0;i<nerrors;i++)y^=1<<epos[i];
/*If there were too many errors, we may not find enough roots to reduce the
syndrome to zero.
We could recompute it to check, but it's much faster just to check that
we have a valid codeword.*/
if(bch15_5_encode(y>>10)==y){
/*Decoding succeeded.*/
*_y=y;
return nerrors;
}
}
/*Decoding failed due to too many bit errors.*/
return -1;
}
unsigned bch15_5_encode(unsigned _x){
return (-(_x&1)&0x0537)^(-(_x>>1&1)&0x0A6E)^(-(_x>>2&1)&0x11EB)^
(-(_x>>3&1)&0x23D6)^(-(_x>>4&1)&0x429B);
}
#if 0
#include <stdio.h>
static unsigned codes[32];
static int hamming(int _a,int _b){
int d;
int n;
d=_a^_b;
for(n=0;d;n++)d&=d-1;
return n;
}
static int closest(int _y){
int min_i;
int min_d;
int i;
int d;
min_i=0;
min_d=hamming(_y,codes[0]);
for(i=1;i<32;i++){
d=hamming(_y,codes[i]);
if(d<min_d){
min_d=d;
min_i=i;
}
}
return codes[min_i];
}
int main(void){
int i;
/*Print a list of the valid (uncorrupt) codewords.*/
for(i=0;i<32;i++)codes[i]=bch15_5_encode(i);
for(i=0;i<32;i++)printf("0x%04X%s",codes[i],i+1<32?" ":"\n");
/*Try to decode all receivable (possibly corrupt) codewords.*/
for(i=0;i<0x8000;i++){
unsigned y;
unsigned z;
int nerrors;
int j;
y=i;
nerrors=bch15_5_correct(&y);
z=closest(i);
if(nerrors<0){
printf("0x%04X->Failed\n",i);
if(hamming(i,z)<=3)printf("Error: 0x%04X should map to 0x%04X\n",i,z);
}
else{
printf("0x%04X->0x%04X\n",i,y);
if(z!=y)printf("Error: 0x%04X should map to 0x%04X\n",i,z);
}
}
return 0;
}
#endif