1 /* dvdisaster: Additional error correction for optical media.
2 * Copyright (C) 2004-2007 Carsten Gnoerlich.
3 * Project home page: http://www.dvdisaster.com
4 * Email: carsten@dvdisaster.com -or- cgnoerlich@fsfe.org
5 *
6 * The Reed-Solomon error correction draws a lot of inspiration - and even code -
7 * from Phil Karn's excellent Reed-Solomon library: http://www.ka9q.net/code/fec/
8 *
9 * This program is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation; either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software
21 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA,
22 * or direct your browser at http://www.gnu.org.
23 */
24
25 #include "l-ec.h"
26 #include <stdlib.h>
27 #include <string.h>
28
29 #include <retro_miscellaneous.h>
30
31 /***
32 *** Mapping between cd frame and parity vectors
33 ***/
34
35 /*
36 * Mapping of frame bytes to P/Q Vectors
37 */
38
PToByteIndex(int p,int i)39 int PToByteIndex(int p, int i)
40 { return 12 + p + i*86;
41 }
42
ByteIndexToP(int b,int * p,int * i)43 void ByteIndexToP(int b, int *p, int *i)
44 { *p = (b-12)%86;
45 *i = (b-12)/86;
46 }
47
QToByteIndex(int q,int i)48 int QToByteIndex(int q, int i)
49 { int offset = 12 + (q & 1);
50
51 if(i == 43) return 2248+q;
52 if(i == 44) return 2300+q;
53
54 q&=~1;
55 return offset + (q*43 + i*88) % 2236;
56 }
57
ByteIndexToQ(int b,int * q,int * i)58 void ByteIndexToQ(int b, int *q, int *i)
59 { int x,y,offset;
60
61 if(b >= 2300)
62 { *i = 44;
63 *q = (b-2300);
64 return;
65 }
66
67 if(b >= 2248)
68 { *i = 43;
69 *q = (b-2248);
70 return;
71 }
72
73 offset = b&1;
74 b = (b-12)/2;
75 x = b/43;
76 y = (b-(x*43))%26;
77 *i = b-(x*43);
78 *q = 2*((x+26-y)%26)+offset;
79 }
80
81 /*
82 * There are 86 vectors of P-parity, yielding a RS(26,24) code.
83 */
84
GetPVector(unsigned char * frame,unsigned char * data,int n)85 void GetPVector(unsigned char *frame, unsigned char *data, int n)
86 { int i;
87 int w_idx = n+12;
88
89 for(i=0; i<26; i++, w_idx+=86)
90 data[i] = frame[w_idx];
91 }
92
SetPVector(unsigned char * frame,unsigned char * data,int n)93 void SetPVector(unsigned char *frame, unsigned char *data, int n)
94 { int i;
95 int w_idx = n+12;
96
97 for(i=0; i<26; i++, w_idx+=86)
98 frame[w_idx] = data[i];
99 }
100
FillPVector(unsigned char * frame,unsigned char data,int n)101 void FillPVector(unsigned char *frame, unsigned char data, int n)
102 { int i;
103 int w_idx = n+12;
104
105 for(i=0; i<26; i++, w_idx+=86)
106 frame[w_idx] = data;
107 }
108
OrPVector(unsigned char * frame,unsigned char value,int n)109 void OrPVector(unsigned char *frame, unsigned char value, int n)
110 { int i;
111 int w_idx = n+12;
112
113 for(i=0; i<26; i++, w_idx+=86)
114 frame[w_idx] |= value;
115 }
116
AndPVector(unsigned char * frame,unsigned char value,int n)117 void AndPVector(unsigned char *frame, unsigned char value, int n)
118 { int i;
119 int w_idx = n+12;
120
121 for(i=0; i<26; i++, w_idx+=86)
122 frame[w_idx] &= value;
123 }
124
125 /*
126 * There are 52 vectors of Q-parity, yielding a RS(45,43) code.
127 */
128
GetQVector(unsigned char * frame,unsigned char * data,int n)129 void GetQVector(unsigned char *frame, unsigned char *data, int n)
130 { int offset = 12 + (n & 1);
131 int w_idx = (n&~1) * 43;
132 int i;
133
134 for(i=0; i<43; i++, w_idx+=88)
135 data[i] = frame[(w_idx % 2236) + offset];
136
137 data[43] = frame[2248 + n];
138 data[44] = frame[2300 + n];
139 }
140
SetQVector(unsigned char * frame,unsigned char * data,int n)141 void SetQVector(unsigned char *frame, unsigned char *data, int n)
142 { int offset = 12 + (n & 1);
143 int w_idx = (n&~1) * 43;
144 int i;
145
146 for(i=0; i<43; i++, w_idx+=88)
147 frame[(w_idx % 2236) + offset] = data[i];
148
149 frame[2248 + n] = data[43];
150 frame[2300 + n] = data[44];
151 }
152
FillQVector(unsigned char * frame,unsigned char data,int n)153 void FillQVector(unsigned char *frame, unsigned char data, int n)
154 { int offset = 12 + (n & 1);
155 int w_idx = (n&~1) * 43;
156 int i;
157
158 for(i=0; i<43; i++, w_idx+=88)
159 frame[(w_idx % 2236) + offset] = data;
160
161 frame[2248 + n] = data;
162 frame[2300 + n] = data;
163 }
164
OrQVector(unsigned char * frame,unsigned char data,int n)165 void OrQVector(unsigned char *frame, unsigned char data, int n)
166 { int offset = 12 + (n & 1);
167 int w_idx = (n&~1) * 43;
168 int i;
169
170 for(i=0; i<43; i++, w_idx+=88)
171 frame[(w_idx % 2236) + offset] |= data;
172
173 frame[2248 + n] |= data;
174 frame[2300 + n] |= data;
175 }
176
AndQVector(unsigned char * frame,unsigned char data,int n)177 void AndQVector(unsigned char *frame, unsigned char data, int n)
178 { int offset = 12 + (n & 1);
179 int w_idx = (n&~1) * 43;
180 int i;
181
182 for(i=0; i<43; i++, w_idx+=88)
183 frame[(w_idx % 2236) + offset] &= data;
184
185 frame[2248 + n] &= data;
186 frame[2300 + n] &= data;
187 }
188
189 /***
190 *** C2 error counting
191 ***/
192
CountC2Errors(unsigned char * frame)193 int CountC2Errors(unsigned char *frame)
194 { int i,count = 0;
195 frame += 2352;
196
197 for(i=0; i<294; i++, frame++)
198 { if(*frame & 0x01) count++;
199 if(*frame & 0x02) count++;
200 if(*frame & 0x04) count++;
201 if(*frame & 0x08) count++;
202 if(*frame & 0x10) count++;
203 if(*frame & 0x20) count++;
204 if(*frame & 0x40) count++;
205 if(*frame & 0x80) count++;
206 }
207
208 return count;
209 }
210
211 /***
212 *** L-EC error correction for CD raw data sectors
213 ***/
214
215 /*
216 * These could be used from ReedSolomonTables,
217 * but hardcoding them is faster.
218 */
219
220 #define NROOTS 2
221 #define LEC_FIRST_ROOT 0 //GF_ALPHA0
222 #define LEC_PRIM_ELEM 1
223 #define LEC_PRIMTH_ROOT 1
224
225 /*
226 * Calculate the error syndrome
227 */
228
DecodePQ(ReedSolomonTables * rt,unsigned char * data,int padding,int * erasure_list,int erasure_count)229 int DecodePQ(ReedSolomonTables *rt, unsigned char *data, int padding,
230 int *erasure_list, int erasure_count)
231 {
232 GaloisTables *gt = rt->gfTables;
233 int syndrome[NROOTS];
234 int lambda[NROOTS+1];
235 int omega[NROOTS+1];
236 int b[NROOTS+1];
237 int reg[NROOTS+1];
238 int root[NROOTS];
239 int loc[NROOTS];
240 int syn_error;
241 int deg_lambda,lambda_roots;
242 int deg_omega;
243 int shortened_size = GF_FIELDMAX - padding;
244 int corrected = 0;
245 int i,j,k;
246 int r,el;
247
248 /*** Form the syndromes: Evaluate data(x) at roots of g(x) */
249
250 for(i=0; i<NROOTS; i++)
251 syndrome[i] = data[0];
252
253 for(j=1; j<shortened_size; j++)
254 for(i=0; i<NROOTS; i++)
255 if(syndrome[i] == 0)
256 syndrome[i] = data[j];
257 else syndrome[i] = data[j] ^ gt->alphaTo[mod_fieldmax(gt->indexOf[syndrome[i]]
258 + (LEC_FIRST_ROOT+i)*LEC_PRIM_ELEM)];
259
260 /*** Convert syndrome to index form, check for nonzero condition. */
261
262 syn_error = 0;
263 for(i=0; i<NROOTS; i++)
264 { syn_error |= syndrome[i];
265 syndrome[i] = gt->indexOf[syndrome[i]];
266 }
267
268 /*** If the syndrome is zero, everything is fine. */
269
270 if(!syn_error)
271 return 0;
272
273 /*** Initialize lambda to be the erasure locator polynomial */
274
275 lambda[0] = 1;
276 lambda[1] = lambda[2] = 0;
277
278 erasure_list[0] += padding;
279 erasure_list[1] += padding;
280
281 if(erasure_count > 2) /* sanity check */
282 erasure_count = 0;
283
284 if(erasure_count > 0)
285 { lambda[1] = gt->alphaTo[mod_fieldmax(LEC_PRIM_ELEM*(GF_FIELDMAX-1-erasure_list[0]))];
286
287 for(i=1; i<erasure_count; i++)
288 { int u = mod_fieldmax(LEC_PRIM_ELEM*(GF_FIELDMAX-1-erasure_list[i]));
289
290 for(j=i+1; j>0; j--)
291 { int tmp = gt->indexOf[lambda[j-1]];
292
293 if(tmp != GF_ALPHA0)
294 lambda[j] ^= gt->alphaTo[mod_fieldmax(u + tmp)];
295 }
296 }
297 }
298
299 for(i=0; i<NROOTS+1; i++)
300 b[i] = gt->indexOf[lambda[i]];
301
302 /*** Berlekamp-Massey algorithm to determine error+erasure locator polynomial */
303
304 r = erasure_count; /* r is the step number */
305 el = erasure_count;
306
307 /* Compute discrepancy at the r-th step in poly-form */
308
309 while(++r <= NROOTS)
310 { int discr_r = 0;
311
312 for(i=0; i<r; i++)
313 if((lambda[i] != 0) && (syndrome[r-i-1] != GF_ALPHA0))
314 discr_r ^= gt->alphaTo[mod_fieldmax(gt->indexOf[lambda[i]] + syndrome[r-i-1])];
315
316 discr_r = gt->indexOf[discr_r];
317
318 if(discr_r == GF_ALPHA0)
319 { /* B(x) = x*B(x) */
320 memmove(b+1, b, NROOTS*sizeof(b[0]));
321 b[0] = GF_ALPHA0;
322 }
323 else
324 { int t[NROOTS+1];
325
326 /* T(x) = lambda(x) - discr_r*x*b(x) */
327 t[0] = lambda[0];
328 for(i=0; i<NROOTS; i++)
329 { if(b[i] != GF_ALPHA0)
330 t[i+1] = lambda[i+1] ^ gt->alphaTo[mod_fieldmax(discr_r + b[i])];
331 else t[i+1] = lambda[i+1];
332 }
333
334 if(2*el <= r+erasure_count-1)
335 { el = r + erasure_count - el;
336
337 /* B(x) <-- inv(discr_r) * lambda(x) */
338 for(i=0; i<=NROOTS; i++)
339 b[i] = (lambda[i] == 0) ? GF_ALPHA0
340 : mod_fieldmax(gt->indexOf[lambda[i]] - discr_r + GF_FIELDMAX);
341 }
342 else
343 { /* 2 lines below: B(x) <-- x*B(x) */
344 memmove(b+1, b, NROOTS*sizeof(b[0]));
345 b[0] = GF_ALPHA0;
346 }
347
348 memcpy(lambda, t, (NROOTS+1)*sizeof(t[0]));
349 }
350 }
351
352 /*** Convert lambda to index form and compute deg(lambda(x)) */
353
354 deg_lambda = 0;
355 for(i=0; i<NROOTS+1; i++)
356 { lambda[i] = gt->indexOf[lambda[i]];
357 if(lambda[i] != GF_ALPHA0)
358 deg_lambda = i;
359 }
360
361 /*** Find roots of the error+erasure locator polynomial by Chien search */
362
363 memcpy(reg+1, lambda+1, NROOTS*sizeof(reg[0]));
364 lambda_roots = 0; /* Number of roots of lambda(x) */
365
366 for(i=1, k=LEC_PRIMTH_ROOT-1; i<=GF_FIELDMAX; i++, k=mod_fieldmax(k+LEC_PRIMTH_ROOT))
367 { int q=1; /* lambda[0] is always 0 */
368
369 for(j=deg_lambda; j>0; j--)
370 { if(reg[j] != GF_ALPHA0)
371 { reg[j] = mod_fieldmax(reg[j] + j);
372 q ^= gt->alphaTo[reg[j]];
373 }
374 }
375
376 if(q != 0) continue; /* Not a root */
377
378 /* store root in index-form and the error location number */
379
380 root[lambda_roots] = i;
381 loc[lambda_roots] = k;
382
383 /* If we've already found max possible roots, abort the search to save time */
384
385 if(++lambda_roots == deg_lambda) break;
386 }
387
388 /* deg(lambda) unequal to number of roots => uncorrectable error detected
389 This is not reliable for very small numbers of roots, e.g. nroots = 2 */
390
391 if(deg_lambda != lambda_roots)
392 { return -1;
393 }
394
395 /* Compute err+eras evaluator poly omega(x) = syn(x)*lambda(x)
396 (modulo x**nroots). in index form. Also find deg(omega). */
397
398 deg_omega = deg_lambda-1;
399
400 for(i=0; i<=deg_omega; i++)
401 { int tmp = 0;
402
403 for(j=i; j>=0; j--)
404 { if((syndrome[i - j] != GF_ALPHA0) && (lambda[j] != GF_ALPHA0))
405 tmp ^= gt->alphaTo[mod_fieldmax(syndrome[i - j] + lambda[j])];
406 }
407
408 omega[i] = gt->indexOf[tmp];
409 }
410
411 /* Compute error values in poly-form.
412 num1 = omega(inv(X(l))),
413 num2 = inv(X(l))**(FIRST_ROOT-1) and
414 den = lambda_pr(inv(X(l))) all in poly-form. */
415
416 for(j=lambda_roots-1; j>=0; j--)
417 { int num1 = 0;
418 int num2;
419 int den;
420 int location = loc[j];
421
422 for(i=deg_omega; i>=0; i--)
423 { if(omega[i] != GF_ALPHA0)
424 num1 ^= gt->alphaTo[mod_fieldmax(omega[i] + i * root[j])];
425 }
426
427 num2 = gt->alphaTo[mod_fieldmax(root[j] * (LEC_FIRST_ROOT - 1) + GF_FIELDMAX)];
428 den = 0;
429
430 /* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
431
432 for(i=MIN(deg_lambda, NROOTS-1) & ~1; i>=0; i-=2)
433 { if(lambda[i+1] != GF_ALPHA0)
434 den ^= gt->alphaTo[mod_fieldmax(lambda[i+1] + i * root[j])];
435 }
436
437 /* Apply error to data */
438
439 if(num1 != 0 && location >= padding)
440 {
441 corrected++;
442 data[location-padding] ^= gt->alphaTo[mod_fieldmax(gt->indexOf[num1] + gt->indexOf[num2]
443 + GF_FIELDMAX - gt->indexOf[den])];
444
445 /* If no erasures were given, at most one error was corrected.
446 Return its position in erasure_list[0]. */
447
448 if(!erasure_count)
449 erasure_list[0] = location-padding;
450 }
451 else
452 return -3;
453 }
454
455 /*** Form the syndromes: Evaluate data(x) at roots of g(x) */
456
457 for(i=0; i<NROOTS; i++)
458 syndrome[i] = data[0];
459
460 for(j=1; j<shortened_size; j++)
461 for(i=0; i<NROOTS; i++)
462 { if(syndrome[i] == 0)
463 syndrome[i] = data[j];
464 else syndrome[i] = data[j] ^ gt->alphaTo[mod_fieldmax(gt->indexOf[syndrome[i]]
465 + (LEC_FIRST_ROOT+i)*LEC_PRIM_ELEM)];
466 }
467
468 /*** Convert syndrome to index form, check for nonzero condition. */
469 for(i=0; i<NROOTS; i++)
470 if(syndrome[i])
471 return -2;
472
473 return corrected;
474 }
475