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 "dvdisaster.h"
26
27 #include "galois-inlines.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 { GaloisTables *gt = rt->gfTables;
232 int syndrome[NROOTS];
233 int lambda[NROOTS+1];
234 int omega[NROOTS+1];
235 int b[NROOTS+1];
236 int reg[NROOTS+1];
237 int root[NROOTS];
238 int loc[NROOTS];
239 int syn_error;
240 int deg_lambda,lambda_roots;
241 int deg_omega;
242 int shortened_size = GF_FIELDMAX - padding;
243 int corrected = 0;
244 int i,j,k;
245 int r,el;
246
247 /*** Form the syndromes: Evaluate data(x) at roots of g(x) */
248
249 for(i=0; i<NROOTS; i++)
250 syndrome[i] = data[0];
251
252 for(j=1; j<shortened_size; j++)
253 for(i=0; i<NROOTS; i++)
254 if(syndrome[i] == 0)
255 syndrome[i] = data[j];
256 else syndrome[i] = data[j] ^ gt->alphaTo[mod_fieldmax(gt->indexOf[syndrome[i]]
257 + (LEC_FIRST_ROOT+i)*LEC_PRIM_ELEM)];
258
259 /*** Convert syndrome to index form, check for nonzero condition. */
260
261 syn_error = 0;
262 for(i=0; i<NROOTS; i++)
263 { syn_error |= syndrome[i];
264 syndrome[i] = gt->indexOf[syndrome[i]];
265 }
266
267 /*** If the syndrome is zero, everything is fine. */
268
269 if(!syn_error)
270 return 0;
271
272 /*** Initialize lambda to be the erasure locator polynomial */
273
274 lambda[0] = 1;
275 lambda[1] = lambda[2] = 0;
276
277 erasure_list[0] += padding;
278 erasure_list[1] += padding;
279
280 if(erasure_count > 2) /* sanity check */
281 erasure_count = 0;
282
283 if(erasure_count > 0)
284 { lambda[1] = gt->alphaTo[mod_fieldmax(LEC_PRIM_ELEM*(GF_FIELDMAX-1-erasure_list[0]))];
285
286 for(i=1; i<erasure_count; i++)
287 { int u = mod_fieldmax(LEC_PRIM_ELEM*(GF_FIELDMAX-1-erasure_list[i]));
288
289 for(j=i+1; j>0; j--)
290 { int tmp = gt->indexOf[lambda[j-1]];
291
292 if(tmp != GF_ALPHA0)
293 lambda[j] ^= gt->alphaTo[mod_fieldmax(u + tmp)];
294 }
295 }
296 }
297
298 for(i=0; i<NROOTS+1; i++)
299 b[i] = gt->indexOf[lambda[i]];
300
301 /*** Berlekamp-Massey algorithm to determine error+erasure locator polynomial */
302
303 r = erasure_count; /* r is the step number */
304 el = erasure_count;
305
306 /* Compute discrepancy at the r-th step in poly-form */
307
308 while(++r <= NROOTS)
309 { int discr_r = 0;
310
311 for(i=0; i<r; i++)
312 if((lambda[i] != 0) && (syndrome[r-i-1] != GF_ALPHA0))
313 discr_r ^= gt->alphaTo[mod_fieldmax(gt->indexOf[lambda[i]] + syndrome[r-i-1])];
314
315 discr_r = gt->indexOf[discr_r];
316
317 if(discr_r == GF_ALPHA0)
318 { /* B(x) = x*B(x) */
319 memmove(b+1, b, NROOTS*sizeof(b[0]));
320 b[0] = GF_ALPHA0;
321 }
322 else
323 { int t[NROOTS+1];
324
325 /* T(x) = lambda(x) - discr_r*x*b(x) */
326 t[0] = lambda[0];
327 for(i=0; i<NROOTS; i++)
328 { if(b[i] != GF_ALPHA0)
329 t[i+1] = lambda[i+1] ^ gt->alphaTo[mod_fieldmax(discr_r + b[i])];
330 else t[i+1] = lambda[i+1];
331 }
332
333 if(2*el <= r+erasure_count-1)
334 { el = r + erasure_count - el;
335
336 /* B(x) <-- inv(discr_r) * lambda(x) */
337 for(i=0; i<=NROOTS; i++)
338 b[i] = (lambda[i] == 0) ? GF_ALPHA0
339 : mod_fieldmax(gt->indexOf[lambda[i]] - discr_r + GF_FIELDMAX);
340 }
341 else
342 { /* 2 lines below: B(x) <-- x*B(x) */
343 memmove(b+1, b, NROOTS*sizeof(b[0]));
344 b[0] = GF_ALPHA0;
345 }
346
347 memcpy(lambda, t, (NROOTS+1)*sizeof(t[0]));
348 }
349 }
350
351 /*** Convert lambda to index form and compute deg(lambda(x)) */
352
353 deg_lambda = 0;
354 for(i=0; i<NROOTS+1; i++)
355 { lambda[i] = gt->indexOf[lambda[i]];
356 if(lambda[i] != GF_ALPHA0)
357 deg_lambda = i;
358 }
359
360 /*** Find roots of the error+erasure locator polynomial by Chien search */
361
362 memcpy(reg+1, lambda+1, NROOTS*sizeof(reg[0]));
363 lambda_roots = 0; /* Number of roots of lambda(x) */
364
365 for(i=1, k=LEC_PRIMTH_ROOT-1; i<=GF_FIELDMAX; i++, k=mod_fieldmax(k+LEC_PRIMTH_ROOT))
366 { int q=1; /* lambda[0] is always 0 */
367
368 for(j=deg_lambda; j>0; j--)
369 { if(reg[j] != GF_ALPHA0)
370 { reg[j] = mod_fieldmax(reg[j] + j);
371 q ^= gt->alphaTo[reg[j]];
372 }
373 }
374
375 if(q != 0) continue; /* Not a root */
376
377 /* store root in index-form and the error location number */
378
379 root[lambda_roots] = i;
380 loc[lambda_roots] = k;
381
382 /* If we've already found max possible roots, abort the search to save time */
383
384 if(++lambda_roots == deg_lambda) break;
385 }
386
387 /* deg(lambda) unequal to number of roots => uncorrectable error detected
388 This is not reliable for very small numbers of roots, e.g. nroots = 2 */
389
390 if(deg_lambda != lambda_roots)
391 { return -1;
392 }
393
394 /* Compute err+eras evaluator poly omega(x) = syn(x)*lambda(x)
395 (modulo x**nroots). in index form. Also find deg(omega). */
396
397 deg_omega = deg_lambda-1;
398
399 for(i=0; i<=deg_omega; i++)
400 { int tmp = 0;
401
402 for(j=i; j>=0; j--)
403 { if((syndrome[i - j] != GF_ALPHA0) && (lambda[j] != GF_ALPHA0))
404 tmp ^= gt->alphaTo[mod_fieldmax(syndrome[i - j] + lambda[j])];
405 }
406
407 omega[i] = gt->indexOf[tmp];
408 }
409
410 /* Compute error values in poly-form.
411 num1 = omega(inv(X(l))),
412 num2 = inv(X(l))**(FIRST_ROOT-1) and
413 den = lambda_pr(inv(X(l))) all in poly-form. */
414
415 for(j=lambda_roots-1; j>=0; j--)
416 { int num1 = 0;
417 int num2;
418 int den;
419 int location = loc[j];
420
421 for(i=deg_omega; i>=0; i--)
422 { if(omega[i] != GF_ALPHA0)
423 num1 ^= gt->alphaTo[mod_fieldmax(omega[i] + i * root[j])];
424 }
425
426 num2 = gt->alphaTo[mod_fieldmax(root[j] * (LEC_FIRST_ROOT - 1) + GF_FIELDMAX)];
427 den = 0;
428
429 /* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
430
431 for(i=MIN(deg_lambda, NROOTS-1) & ~1; i>=0; i-=2)
432 { if(lambda[i+1] != GF_ALPHA0)
433 den ^= gt->alphaTo[mod_fieldmax(lambda[i+1] + i * root[j])];
434 }
435
436 /* Apply error to data */
437
438 if(num1 != 0 && location >= padding)
439 {
440 corrected++;
441 data[location-padding] ^= gt->alphaTo[mod_fieldmax(gt->indexOf[num1] + gt->indexOf[num2]
442 + GF_FIELDMAX - gt->indexOf[den])];
443
444 /* If no erasures were given, at most one error was corrected.
445 Return its position in erasure_list[0]. */
446
447 if(!erasure_count)
448 erasure_list[0] = location-padding;
449 }
450 #if 1
451 else return -3;
452 #endif
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 #if 1
470 for(i=0; i<NROOTS; i++)
471 if(syndrome[i])
472 return -2;
473 #endif
474
475 return corrected;
476 }
477
478
479