Line data Source code
1 : /* General purpose Reed-Solomon codec for 8-bit symbols or less
2 : * Copyright 2003 Phil Karn, KA9Q
3 : * May be used under the terms of the GNU Lesser General Public License (LGPL)
4 : *
5 : * This version tweaked by Philip Heron <phil@sanslogic.co.uk>
6 : */
7 :
8 : #include <string.h>
9 : #include "ReedSolomon.h"
10 :
11 : static const uint8_t ALPHA_TO[] = {
12 : 0x01,0x02,0x04,0x08,0x10,0x20,0x40,0x80,0x87,0x89,0x95,0xAD,0xDD,0x3D,0x7A,0xF4,
13 : 0x6F,0xDE,0x3B,0x76,0xEC,0x5F,0xBE,0xFB,0x71,0xE2,0x43,0x86,0x8B,0x91,0xA5,0xCD,
14 : 0x1D,0x3A,0x74,0xE8,0x57,0xAE,0xDB,0x31,0x62,0xC4,0x0F,0x1E,0x3C,0x78,0xF0,0x67,
15 : 0xCE,0x1B,0x36,0x6C,0xD8,0x37,0x6E,0xDC,0x3F,0x7E,0xFC,0x7F,0xFE,0x7B,0xF6,0x6B,
16 : 0xD6,0x2B,0x56,0xAC,0xDF,0x39,0x72,0xE4,0x4F,0x9E,0xBB,0xF1,0x65,0xCA,0x13,0x26,
17 : 0x4C,0x98,0xB7,0xE9,0x55,0xAA,0xD3,0x21,0x42,0x84,0x8F,0x99,0xB5,0xED,0x5D,0xBA,
18 : 0xF3,0x61,0xC2,0x03,0x06,0x0C,0x18,0x30,0x60,0xC0,0x07,0x0E,0x1C,0x38,0x70,0xE0,
19 : 0x47,0x8E,0x9B,0xB1,0xE5,0x4D,0x9A,0xB3,0xE1,0x45,0x8A,0x93,0xA1,0xC5,0x0D,0x1A,
20 : 0x34,0x68,0xD0,0x27,0x4E,0x9C,0xBF,0xF9,0x75,0xEA,0x53,0xA6,0xCB,0x11,0x22,0x44,
21 : 0x88,0x97,0xA9,0xD5,0x2D,0x5A,0xB4,0xEF,0x59,0xB2,0xE3,0x41,0x82,0x83,0x81,0x85,
22 : 0x8D,0x9D,0xBD,0xFD,0x7D,0xFA,0x73,0xE6,0x4B,0x96,0xAB,0xD1,0x25,0x4A,0x94,0xAF,
23 : 0xD9,0x35,0x6A,0xD4,0x2F,0x5E,0xBC,0xFF,0x79,0xF2,0x63,0xC6,0x0B,0x16,0x2C,0x58,
24 : 0xB0,0xE7,0x49,0x92,0xA3,0xC1,0x05,0x0A,0x14,0x28,0x50,0xA0,0xC7,0x09,0x12,0x24,
25 : 0x48,0x90,0xA7,0xC9,0x15,0x2A,0x54,0xA8,0xD7,0x29,0x52,0xA4,0xCF,0x19,0x32,0x64,
26 : 0xC8,0x17,0x2E,0x5C,0xB8,0xF7,0x69,0xD2,0x23,0x46,0x8C,0x9F,0xB9,0xF5,0x6D,0xDA,
27 : 0x33,0x66,0xCC,0x1F,0x3E,0x7C,0xF8,0x77,0xEE,0x5B,0xB6,0xEB,0x51,0xA2,0xC3,0x00,
28 : };
29 :
30 : static const uint8_t INDEX_OF[] = {
31 : 0xFF,0x00,0x01,0x63,0x02,0xC6,0x64,0x6A,0x03,0xCD,0xC7,0xBC,0x65,0x7E,0x6B,0x2A,
32 : 0x04,0x8D,0xCE,0x4E,0xC8,0xD4,0xBD,0xE1,0x66,0xDD,0x7F,0x31,0x6C,0x20,0x2B,0xF3,
33 : 0x05,0x57,0x8E,0xE8,0xCF,0xAC,0x4F,0x83,0xC9,0xD9,0xD5,0x41,0xBE,0x94,0xE2,0xB4,
34 : 0x67,0x27,0xDE,0xF0,0x80,0xB1,0x32,0x35,0x6D,0x45,0x21,0x12,0x2C,0x0D,0xF4,0x38,
35 : 0x06,0x9B,0x58,0x1A,0x8F,0x79,0xE9,0x70,0xD0,0xC2,0xAD,0xA8,0x50,0x75,0x84,0x48,
36 : 0xCA,0xFC,0xDA,0x8A,0xD6,0x54,0x42,0x24,0xBF,0x98,0x95,0xF9,0xE3,0x5E,0xB5,0x15,
37 : 0x68,0x61,0x28,0xBA,0xDF,0x4C,0xF1,0x2F,0x81,0xE6,0xB2,0x3F,0x33,0xEE,0x36,0x10,
38 : 0x6E,0x18,0x46,0xA6,0x22,0x88,0x13,0xF7,0x2D,0xB8,0x0E,0x3D,0xF5,0xA4,0x39,0x3B,
39 : 0x07,0x9E,0x9C,0x9D,0x59,0x9F,0x1B,0x08,0x90,0x09,0x7A,0x1C,0xEA,0xA0,0x71,0x5A,
40 : 0xD1,0x1D,0xC3,0x7B,0xAE,0x0A,0xA9,0x91,0x51,0x5B,0x76,0x72,0x85,0xA1,0x49,0xEB,
41 : 0xCB,0x7C,0xFD,0xC4,0xDB,0x1E,0x8B,0xD2,0xD7,0x92,0x55,0xAA,0x43,0x0B,0x25,0xAF,
42 : 0xC0,0x73,0x99,0x77,0x96,0x5C,0xFA,0x52,0xE4,0xEC,0x5F,0x4A,0xB6,0xA2,0x16,0x86,
43 : 0x69,0xC5,0x62,0xFE,0x29,0x7D,0xBB,0xCC,0xE0,0xD3,0x4D,0x8C,0xF2,0x1F,0x30,0xDC,
44 : 0x82,0xAB,0xE7,0x56,0xB3,0x93,0x40,0xD8,0x34,0xB0,0xEF,0x26,0x37,0x0C,0x11,0x44,
45 : 0x6F,0x78,0x19,0x9A,0x47,0x74,0xA7,0xC1,0x23,0x53,0x89,0xFB,0x14,0x5D,0xF8,0x97,
46 : 0x2E,0x4B,0xB9,0x60,0x0F,0xED,0x3E,0xE5,0xF6,0x87,0xA5,0x17,0x3A,0xA3,0x3C,0xB7,
47 : };
48 :
49 : static const uint8_t GENPOLY[] = {
50 : 0x00,0xF9,0x3B,0x42,0x04,0x2B,0x7E,0xFB,0x61,0x1E,0x03,0xD5,0x32,0x42,0xAA,0x05,
51 : 0x18,0x05,0xAA,0x42,0x32,0xD5,0x03,0x1E,0x61,0xFB,0x7E,0x2B,0x04,0x42,0x3B,0xF9,
52 : 0x00,
53 : };
54 :
55 75103 : static inline int mod255(int x)
56 : {
57 132382 : while(x >= 255)
58 : {
59 57279 : x -= 255;
60 57279 : x = (x >> 8) + (x & 255);
61 : }
62 75103 : return(x);
63 : }
64 : #define MODNN(x) mod255(x)
65 :
66 : #define MM (8)
67 : #define NN (255)
68 : #define NROOTS (32)
69 : #define FCR (112)
70 : #define PRIM (11)
71 : #define IPRIM (116)
72 :
73 : #define MIN(a,b) ((a) < (b) ? (a) : (b))
74 : #define A0 (NN) /* Special reserved value encoding zero in index form */
75 :
76 : /* Portable C version */
77 3 : void rlb_encode_rs_8(const uint8_t *data, uint8_t *parity, int pad)
78 : {
79 : int i, j;
80 :
81 3 : memset(parity, 0, NROOTS * sizeof(uint8_t));
82 :
83 662 : for(i = 0; i < NN - NROOTS - pad; i++)
84 : {
85 659 : uint8_t feedback = INDEX_OF[data[i] ^ parity[0]];
86 659 : if(feedback != A0) /* feedback term is non-zero */
87 : {
88 21088 : for(j = 1; j < NROOTS; j++)
89 20429 : parity[j] ^= ALPHA_TO[mod255(feedback + GENPOLY[NROOTS - j])];
90 : }
91 :
92 : /* Shift */
93 659 : memmove(&parity[0], &parity[1], sizeof(uint8_t) * (NROOTS - 1));
94 659 : if(feedback != A0)
95 659 : parity[NROOTS - 1] = ALPHA_TO[mod255(feedback + GENPOLY[0])];
96 : else
97 0 : parity[NROOTS - 1] = 0;
98 : }
99 3 : }
100 :
101 6 : int rlb_decode_rs_8(uint8_t *data, int *eras_pos, int no_eras, int pad)
102 : {
103 : int deg_lambda, el, deg_omega;
104 : int i, j, r, k;
105 : uint8_t u, q, tmp, num1, num2, den, discr_r;
106 : uint8_t lambda[NROOTS + 1], s[NROOTS]; /* Err+Eras Locator poly
107 : * and syndrome poly */
108 : uint8_t b[NROOTS + 1], t[NROOTS + 1], omega[NROOTS + 1];
109 : uint8_t root[NROOTS], reg[NROOTS + 1], loc[NROOTS];
110 : int syn_error, count;
111 :
112 6 : if(pad < 0 || pad > 222) return(-1);
113 :
114 : /* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
115 198 : for(i = 0; i < NROOTS; i++) s[i] = data[0];
116 :
117 1520 : for(j = 1; j < NN - pad; j++)
118 : {
119 49962 : for(i = 0; i < NROOTS; i++)
120 : {
121 48448 : if(s[i] == 0) s[i] = data[j];
122 48161 : else s[i] = data[j] ^ ALPHA_TO[MODNN(INDEX_OF[s[i]] + (FCR + i) * PRIM)];
123 : }
124 : }
125 :
126 : /* Convert syndromes to index form, checking for nonzero condition */
127 6 : syn_error = 0;
128 198 : for(i = 0; i < NROOTS; i++)
129 : {
130 192 : syn_error |= s[i];
131 192 : s[i] = INDEX_OF[s[i]];
132 : }
133 :
134 6 : if(!syn_error)
135 : {
136 : /* if syndrome is zero, data[] is a codeword and there are no
137 : * errors to correct. So return data[] unmodified
138 : */
139 2 : count = 0;
140 2 : goto finish;
141 : }
142 :
143 4 : memset(&lambda[1], 0, NROOTS * sizeof(lambda[0]));
144 4 : lambda[0] = 1;
145 :
146 4 : if(no_eras > 0)
147 : {
148 : /* Init lambda to be the erasure locator polynomial */
149 1 : lambda[1] = ALPHA_TO[MODNN(PRIM * (NN - 1 - eras_pos[0]))];
150 1 : for(i = 1; i < no_eras; i++)
151 : {
152 0 : u = MODNN(PRIM * (NN - 1 - eras_pos[i]));
153 0 : for(j = i + 1; j > 0; j--)
154 : {
155 0 : tmp = INDEX_OF[lambda[j - 1]];
156 0 : if(tmp != A0) lambda[j] ^= ALPHA_TO[MODNN(u + tmp)];
157 : }
158 : }
159 :
160 : }
161 :
162 136 : for(i = 0; i < NROOTS + 1; i++)
163 132 : b[i] = INDEX_OF[lambda[i]];
164 :
165 : /*
166 : * Begin Berlekamp-Massey algorithm to determine error+erasure
167 : * locator polynomial
168 : */
169 4 : r = no_eras;
170 4 : el = no_eras;
171 131 : while(++r <= NROOTS) /* r is the step number */
172 : {
173 : /* Compute discrepancy at the r-th step in poly-form */
174 127 : discr_r = 0;
175 2238 : for(i = 0; i < r; i++)
176 : {
177 2111 : if((lambda[i] != 0) && (s[r - i - 1] != A0))
178 : {
179 456 : discr_r ^= ALPHA_TO[MODNN(INDEX_OF[lambda[i]] + s[r - i - 1])];
180 : }
181 : }
182 127 : discr_r = INDEX_OF[discr_r]; /* Index form */
183 127 : if(discr_r == A0)
184 : {
185 : /* 2 lines below: B(x) <-- x*B(x) */
186 92 : memmove(&b[1], b, NROOTS * sizeof(b[0]));
187 92 : b[0] = A0;
188 : }
189 : else
190 : {
191 : /* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
192 35 : t[0] = lambda[0];
193 1155 : for(i = 0; i < NROOTS; i++)
194 : {
195 1120 : if(b[i] != A0)
196 248 : t[i + 1] = lambda[i + 1] ^ ALPHA_TO[MODNN(discr_r + b[i])];
197 : else
198 872 : t[i + 1] = lambda[i + 1];
199 : }
200 :
201 35 : if(2 * el <= r + no_eras - 1)
202 : {
203 17 : el = r + no_eras - el;
204 : /*
205 : * 2 lines below: B(x) <-- inv(discr_r) *
206 : * lambda(x)
207 : */
208 578 : for(i = 0; i <= NROOTS; i++)
209 561 : b[i] = (lambda[i] == 0) ? A0 : MODNN(INDEX_OF[lambda[i]] - discr_r + NN);
210 : }
211 : else
212 : {
213 : /* 2 lines below: B(x) <-- x*B(x) */
214 18 : memmove(&b[1], b, NROOTS * sizeof(b[0]));
215 18 : b[0] = A0;
216 : }
217 :
218 35 : memcpy(lambda, t, (NROOTS + 1) * sizeof(t[0]));
219 : }
220 : }
221 :
222 : /* Convert lambda to index form and compute deg(lambda(x)) */
223 4 : deg_lambda = 0;
224 136 : for(i = 0; i < NROOTS + 1; i++)
225 : {
226 132 : lambda[i] = INDEX_OF[lambda[i]];
227 132 : if(lambda[i] != A0) deg_lambda = i;
228 : }
229 :
230 : /* Find roots of the error+erasure locator polynomial by Chien search */
231 4 : memcpy(®[1], &lambda[1], NROOTS * sizeof(reg[0]));
232 4 : count = 0; /* Number of roots of lambda(x) */
233 594 : for(i = 1, k = IPRIM - 1; i <= NN; i++, k = MODNN(k + IPRIM))
234 : {
235 593 : q = 1; /* lambda[0] is always 0 */
236 5011 : for(j = deg_lambda; j > 0; j--)
237 : {
238 4418 : if(reg[j] != A0)
239 : {
240 4418 : reg[j] = MODNN(reg[j] + j);
241 4418 : q ^= ALPHA_TO[reg[j]];
242 : }
243 : }
244 :
245 593 : if (q != 0) continue; /* Not a root */
246 :
247 : /* store root (index-form) and error location number */
248 5 : root[count] = i;
249 5 : loc[count] = k;
250 : /* If we've already found max possible roots,
251 : * abort the search to save time
252 : */
253 5 : if(++count == deg_lambda) break;
254 : }
255 :
256 4 : if(deg_lambda != count)
257 : {
258 : /*
259 : * deg(lambda) unequal to number of roots => uncorrectable
260 : * error detected
261 : */
262 1 : count = -1;
263 1 : goto finish;
264 : }
265 :
266 : /*
267 : * Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
268 : * x**NROOTS). in index form. Also find deg(omega).
269 : */
270 3 : deg_omega = deg_lambda - 1;
271 6 : for(i = 0; i <= deg_omega; i++)
272 : {
273 3 : tmp = 0;
274 6 : for(j = i; j >= 0; j--)
275 : {
276 3 : if((s[i - j] != A0) && (lambda[j] != A0))
277 3 : tmp ^= ALPHA_TO[MODNN(s[i - j] + lambda[j])];
278 : }
279 3 : omega[i] = INDEX_OF[tmp];
280 : }
281 :
282 : /*
283 : * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
284 : * inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form
285 : */
286 6 : for(j = count - 1; j >= 0; j--)
287 : {
288 3 : num1 = 0;
289 6 : for(i = deg_omega; i >= 0; i--)
290 : {
291 3 : if(omega[i] != A0) num1 ^= ALPHA_TO[MODNN(omega[i] + i * root[j])];
292 : }
293 3 : num2 = ALPHA_TO[MODNN(root[j] * (FCR - 1) + NN)];
294 3 : den = 0;
295 :
296 : /* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
297 6 : for(i = MIN(deg_lambda, NROOTS - 1) & ~1; i >= 0; i -= 2)
298 : {
299 3 : if(lambda[i + 1] != A0) den ^= ALPHA_TO[MODNN(lambda[i + 1] + i * root[j])];
300 : }
301 :
302 : /* Apply error to data */
303 3 : if(num1 != 0 && loc[j] >= pad)
304 : {
305 3 : data[loc[j] - pad] ^= ALPHA_TO[MODNN(INDEX_OF[num1] + INDEX_OF[num2] + NN - INDEX_OF[den])];
306 : }
307 : }
308 :
309 3 : finish:
310 6 : if(eras_pos != NULL)
311 : {
312 2 : for(i = 0; i < count; i++) eras_pos[i] = loc[i];
313 : }
314 :
315 6 : return(count);
316 : }
|
