кусок библиотеки восстановления ошибок с помощью кодов Рида-Соломона

 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
// полагаем, что нам передали паекты с кол-вом ошибок <= (n-k)
// иначе... падаем по переполнению стэка ;)
void rs_decode(u_int8_t *pkt[RS_PKT], u_int8_t bugmask[RS_PKT], int pkt_size, int k) {
int errndx[RS_PKT-k]; // больше мы не восстановим
int errnum;
int slice, i, j;
errnum = 0;
for (i = 0; i < RS_PKT; i++) {
if (bugmask[i]) {
memset(pkt[i], 0, pkt_size); // GF_NUL
errndx[errnum] = i;
errnum++;
}
}
// индексы ошибок заполнены
for (slice = 0; slice < pkt_size; slice++) {
u_int8_t a[errnum][errnum+1]; // первый индекс - номер строки, второй - столбца
u_int8_t r[RS_PKT]; // принятый полином, нарезанный из слайса
for (i = 0; i < RS_PKT; i++)
r[i] = pkt[i][slice];
for (i = 0; i < errnum; i++) { // заполняем систему уравнений
for (j = 0; j < errnum; j++) { // заполняем основную матрицу системы
a[i][j] = gf_pow(i+1, errndx[j]);
}
a[i][errnum] = gf_poly_subst(r, RS_PKT-1, i+1); // и столбец свободных членов
}
// приводим матрицу к диагональному виду методом Гаусса
for (i = 0; i < errnum; i++) { // к верхнетреугольному, i - строка
int line;
for (line = i; line < errnum; line++) { // найдем строчку с ненулевым коэфф-ом
if (a[line][i] != GF_NUL)
break;
}
if (line != i) { // i-я строчка содержит нуль на i-м месте
u_int8_t foo[errnum+1]; // temporary storage for line
memcpy(foo, a[line], errnum+1);
memcpy(a[line], a[i], errnum+1);
memcpy(a[i], foo, errnum+1);
}
if (a[i][i] != GF_MUL_ONE) { // если первый член != 1, приводим его
u_int8_t divisor = a[i][i];
for (j = i; j < errnum+1; j++) {
a[i][j] = gf_div(a[i][j], divisor);
}
}
for (line = i+1; line < errnum; line++) { // вычитаем из низлежащих i-ю строку
u_int8_t coef = a[line][i];
for (j = i; j < errnum+1; j++) {
a[line][j] = gf_add( a[line][j], gf_mul(a[i][j], coef) );
}
}
}
for (i = errnum-1; i >= 0; i--) { // к диагональному виду, i - текущая строка
for (j = 0; j < i; j++) { // вычесть её из каждой строчки
a[j][errnum] = gf_add( a[j][errnum], gf_mul( a[i][errnum], a[j][i]) );
a[j][i] = GF_NUL;
}
}
//~ for (i = 0; i < errnum; i++) {
//~ for (j = 0; j < errnum+1; j++) {
//~ printf("%3.3i ", a[i][j]);
//~ }
//~ printf("\n");
//~ }
//~ #warning test test test
for (i = 0; i < errnum; i++)
pkt[errndx[i]][slice] = a[i][i];
}
}