1 // This file is part of par2cmdline (a PAR 2.0 compatible file verification and
2 // repair tool). See http://parchive.sourceforge.net for details of PAR 2.0.
3 //
4 // Copyright (c) 2003 Peter Brian Clements
5 // Copyright (c) 2019 Michael D. Nahas
6 //
7 // par2cmdline is free software; you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation; either version 2 of the License, or
10 // (at your option) any later version.
11 //
12 // par2cmdline is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with this program; if not, write to the Free Software
19 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20
21 #ifndef __REEDSOLOMON_H__
22 #define __REEDSOLOMON_H__
23
24 // The ReedSolomon object is used to calculate and store the matrix
25 // used during recovery block creation or data block reconstruction.
26 //
27 // During initialisation, one RSOutputRow object is created for each
28 // recovery block that either needs to be created or is available for
29 // use.
30
31 class RSOutputRow
32 {
33 public:
RSOutputRow(void)34 RSOutputRow(void) {};
RSOutputRow(bool _present,u16 _exponent)35 RSOutputRow(bool _present, u16 _exponent) : present(_present), exponent(_exponent) {}
36
37 public:
38 bool present;
39 u16 exponent;
40 };
41
42 template<class g>
43 class ReedSolomon
44 {
45 public:
46 typedef g G;
47
48 ReedSolomon(void);
49 ~ReedSolomon(void);
50
51 // Set which input blocks are present or missing
52 // Some input blocks are present
53 bool SetInput(const vector<bool> &present, std::ostream &sout, std::ostream &serr);
54 // All input blocks are present
55 bool SetInput(u32 count, std::ostream &sout, std::ostream &serr);
56
57 // Set which output block are available or need to be computed
58 bool SetOutput(bool present, u16 exponent);
59 bool SetOutput(bool present, u16 lowexponent, u16 highexponent);
60
61 // Compute the RS Matrix
62 bool Compute(NoiseLevel noiselevel, std::ostream &sout, std::ostream &serr);
63
64 // Process a block of data
65 bool Process(size_t size, // The size of the block of data
66 u32 inputindex, // The column in the RS matrix
67 const void *inputbuffer, // Buffer containing input data
68 u32 outputindex, // The row in the RS matrix
69 void *outputbuffer); // Buffer containing output data
70 private:
71 bool InternalProcess(const g &factor, size_t size, const void *inputbuffer, void *outputbuffer); // Optimization
72
73 protected:
74 // Perform Gaussian Elimination
75 bool GaussElim(NoiseLevel noiselevel,
76 std::ostream &sout,
77 std::ostream &serr,
78 unsigned int rows,
79 unsigned int leftcols,
80 G *leftmatrix,
81 G *rightmatrix,
82 unsigned int datamissing);
83
84 protected:
85 u32 inputcount; // Total number of input blocks
86
87 u32 datapresent; // Number of input blocks that are present
88 u32 datamissing; // Number of input blocks that are missing
89 u32 *datapresentindex; // The index numbers of the data blocks that are present
90 u32 *datamissingindex; // The index numbers of the data blocks that are missing
91
92 typename G::ValueType *database;// The "base" value to use for each input block
93
94 u32 outputcount; // Total number of output blocks
95
96 u32 parpresent; // Number of output blocks that are present
97 u32 parmissing; // Number of output blocks that are missing
98 u32 *parpresentindex; // The index numbers of the output blocks that are present
99 u32 *parmissingindex; // The index numbers of the output blocks that are missing
100
101 vector<RSOutputRow> outputrows; // Details of the output blocks
102
103 G *leftmatrix; // The main matrix
104
105 // When the matrices are initialised: values of the form base ^ exponent are
106 // stored (where the base values are obtained from database[] and the exponent
107 // values are obtained from outputrows[]).
108
109 #ifdef LONGMULTIPLY
110 GaloisLongMultiplyTable<g> *glmt; // A multiplication table used by Process()
111 #endif
112
113 private:
114 // private copy constructor to prevent any misuse.
115 ReedSolomon(const ReedSolomon &);
116 ReedSolomon& operator=(const ReedSolomon &);
117 };
118
119 template<class g>
ReedSolomon(void)120 inline ReedSolomon<g>::ReedSolomon(void)
121 : outputrows()
122 {
123 inputcount = 0;
124
125 datapresent = 0;
126 datamissing = 0;
127 datapresentindex = 0;
128 datamissingindex = 0;
129 database = 0;
130
131 outputcount = 0;
132
133 parpresent = 0;
134 parmissing = 0;
135 parpresentindex = 0;
136 parmissingindex = 0;
137
138 leftmatrix = 0;
139
140 #ifdef LONGMULTIPLY
141 glmt = new GaloisLongMultiplyTable<g>;
142 #endif
143 }
144
145 template<class g>
~ReedSolomon(void)146 inline ReedSolomon<g>::~ReedSolomon(void)
147 {
148 delete [] datapresentindex;
149 delete [] datamissingindex;
150 delete [] database;
151 delete [] parpresentindex;
152 delete [] parmissingindex;
153 delete [] leftmatrix;
154
155 #ifdef LONGMULTIPLY
156 delete glmt;
157 #endif
158 }
159
160 template<class g>
Process(size_t size,u32 inputindex,const void * inputbuffer,u32 outputindex,void * outputbuffer)161 inline bool ReedSolomon<g>::Process(size_t size, u32 inputindex, const void *inputbuffer, u32 outputindex, void *outputbuffer)
162 {
163 // Optimization: it occurs frequently the function exits early on, so inline the start.
164 // This resulted in a speed gain of approx. 8% in repairing.
165
166 // Look up the appropriate element in the RS matrix
167 g factor = leftmatrix[outputindex * (datapresent + datamissing) + inputindex];
168 // Do nothing if the factor happens to be 0
169 if (factor == 0)
170 return eSuccess;
171 return this->InternalProcess (factor, size, inputbuffer, outputbuffer);
172 }
173
174 u32 gcd(u32 a, u32 b);
175
176 // Record whether the recovery block with the specified
177 // exponent values is present or missing.
178 template<class g>
SetOutput(bool present,u16 exponent)179 inline bool ReedSolomon<g>::SetOutput(bool present, u16 exponent)
180 {
181 // Store the exponent and whether or not the recovery block is present or missing
182 outputrows.push_back(RSOutputRow(present, exponent));
183
184 outputcount++;
185
186 // Update the counts.
187 if (present)
188 {
189 parpresent++;
190 }
191 else
192 {
193 parmissing++;
194 }
195
196 return true;
197 }
198
199 // Record whether the recovery blocks with the specified
200 // range of exponent values are present or missing.
201 template<class g>
SetOutput(bool present,u16 lowexponent,u16 highexponent)202 inline bool ReedSolomon<g>::SetOutput(bool present, u16 lowexponent, u16 highexponent)
203 {
204 for (unsigned int exponent=lowexponent; exponent<=highexponent; exponent++)
205 {
206 if (!SetOutput(present, exponent))
207 return false;
208 }
209
210 return true;
211 }
212
213 // Construct the Vandermonde matrix and solve it if necessary
214 template<class g>
Compute(NoiseLevel noiselevel,std::ostream & sout,std::ostream & serr)215 inline bool ReedSolomon<g>::Compute(NoiseLevel noiselevel, std::ostream &sout, std::ostream &serr)
216 {
217 u32 outcount = datamissing + parmissing;
218 u32 incount = datapresent + datamissing;
219
220 if (datamissing > parpresent)
221 {
222 serr << "Not enough recovery blocks." << endl;
223 return false;
224 }
225 else if (outcount == 0)
226 {
227 serr << "No output blocks." << endl;
228 return false;
229 }
230
231 if (noiselevel > nlQuiet)
232 sout << "Computing Reed Solomon matrix." << endl;
233
234 /* Layout of RS Matrix:
235 NOTE: The second set of columns represents the parity vectors present,
236 but this only uses datamissing of them. Otherwise, it is over constrained.
237
238 datamissing
239 datapresent <=parpresent datamissing parmissing
240 / | \ / | \
241 | (ppi[row])| | | (ppi[row])| |
242 datamissing | ^ | I | | ^ | 0 |
243 |(dpi[col]) | | |(dmi[col]) | |
244 +---------------------+-------------+ +---------------------+-----------+
245 | (pmi[row])| | | (pmi[row])| |
246 parmissing | ^ | 0 | | ^ | I |
247 |(dpi[col]) | | |(dmi[col]) | |
248 \ | / \ | /
249 */
250
251 // Allocate the left hand matrix
252
253 leftmatrix = new G[outcount * incount];
254 memset(leftmatrix, 0, outcount * incount * sizeof(G));
255
256 // Allocate the right hand matrix only if we are recovering
257
258 G *rightmatrix = 0;
259 if (datamissing > 0)
260 {
261 rightmatrix = new G[outcount * outcount];
262 memset(rightmatrix, 0, outcount *outcount * sizeof(G));
263 }
264
265 // Fill in the two matrices:
266
267 vector<RSOutputRow>::const_iterator outputrow = outputrows.begin();
268
269 // One row for each present recovery block that will be used for a missing data block
270 for (unsigned int row=0; row<datamissing; row++)
271 {
272 // Define MPDL to skip reporting and speed things up
273 #ifndef MPDL
274 if (noiselevel > nlQuiet)
275 {
276 int progress = row * 1000 / (datamissing+parmissing);
277 sout << "Constructing: " << progress/10 << '.' << progress%10 << "%\r" << flush;
278 }
279 #endif
280
281 // Get the exponent of the next present recovery block
282 while (!outputrow->present)
283 {
284 outputrow++;
285 }
286 u16 exponent = outputrow->exponent;
287
288 // One column for each present data block
289 for (unsigned int col=0; col<datapresent; col++)
290 {
291 leftmatrix[row * incount + col] = G(database[datapresentindex[col]]).pow(exponent);
292 }
293 // One column for each each present recovery block that will be used for a missing data block
294 for (unsigned int col=0; col<datamissing; col++)
295 {
296 leftmatrix[row * incount + col + datapresent] = (row == col) ? 1 : 0;
297 }
298
299 if (datamissing > 0)
300 {
301 // One column for each missing data block
302 for (unsigned int col=0; col<datamissing; col++)
303 {
304 rightmatrix[row * outcount + col] = G(database[datamissingindex[col]]).pow(exponent);
305 }
306 // One column for each missing recovery block
307 for (unsigned int col=0; col<parmissing; col++)
308 {
309 rightmatrix[row * outcount + col + datamissing] = 0;
310 }
311 }
312
313 outputrow++;
314 }
315 // One row for each recovery block being computed
316 outputrow = outputrows.begin();
317 for (unsigned int row=0; row<parmissing; row++)
318 {
319 // Define MPDL to skip reporting and speed things up
320 #ifndef MPDL
321 if (noiselevel > nlQuiet)
322 {
323 int progress = (row+datamissing) * 1000 / (datamissing+parmissing);
324 sout << "Constructing: " << progress/10 << '.' << progress%10 << "%\r" << flush;
325 }
326 #endif
327
328 // Get the exponent of the next missing recovery block
329 while (outputrow->present)
330 {
331 outputrow++;
332 }
333 u16 exponent = outputrow->exponent;
334
335 // One column for each present data block
336 for (unsigned int col=0; col<datapresent; col++)
337 {
338 leftmatrix[(row+datamissing) * incount + col] = G(database[datapresentindex[col]]).pow(exponent);
339 }
340 // One column for each each present recovery block that will be used for a missing data block
341 for (unsigned int col=0; col<datamissing; col++)
342 {
343 leftmatrix[(row+datamissing) * incount + col + datapresent] = 0;
344 }
345
346 if (datamissing > 0)
347 {
348 // One column for each missing data block
349 for (unsigned int col=0; col<datamissing; col++)
350 {
351 rightmatrix[(row+datamissing) * outcount + col] = G(database[datamissingindex[col]]).pow(exponent);
352 }
353 // One column for each missing recovery block
354 for (unsigned int col=0; col<parmissing; col++)
355 {
356 rightmatrix[(row+datamissing) * outcount + col + datamissing] = (row == col) ? 1 : 0;
357 }
358 }
359
360 outputrow++;
361 }
362 if (noiselevel > nlQuiet)
363 sout << "Constructing: done." << endl;
364
365 // Solve the matrices only if recovering data
366 if (datamissing > 0)
367 {
368 // Perform Gaussian Elimination and then delete the right matrix (which
369 // will no longer be required).
370 bool success = GaussElim(noiselevel, sout, serr, outcount, incount, leftmatrix, rightmatrix, datamissing);
371 delete [] rightmatrix;
372 return success;
373 }
374
375 return true;
376 }
377
378 // Use Gaussian Elimination to solve the matrices
379 template<class g>
GaussElim(NoiseLevel noiselevel,std::ostream & sout,std::ostream & serr,unsigned int rows,unsigned int leftcols,G * leftmatrix,G * rightmatrix,unsigned int datamissing)380 inline bool ReedSolomon<g>::GaussElim(NoiseLevel noiselevel, std::ostream &sout, std::ostream &serr, unsigned int rows, unsigned int leftcols, G *leftmatrix, G *rightmatrix, unsigned int datamissing)
381 {
382 if (noiselevel >= nlDebug)
383 {
384 for (unsigned int row=0; row<rows; row++)
385 {
386 sout << ((row==0) ? "/" : (row==rows-1) ? "\\" : "|");
387 for (unsigned int col=0; col<leftcols; col++)
388 {
389 sout << " "
390 << hex << setw(G::Bits>8?4:2) << setfill('0')
391 << (unsigned int)leftmatrix[row*leftcols+col];
392 }
393 sout << ((row==0) ? " \\ /" : (row==rows-1) ? " / \\" : " | |");
394 for (unsigned int col=0; col<rows; col++)
395 {
396 sout << " "
397 << hex << setw(G::Bits>8?4:2) << setfill('0')
398 << (unsigned int)rightmatrix[row*rows+col];
399 }
400 sout << ((row==0) ? " \\" : (row==rows-1) ? " /" : " | |");
401 sout << endl;
402
403 sout << dec << setw(0) << setfill(' ');
404 }
405 }
406
407 // Because the matrices being operated on are Vandermonde matrices
408 // they are guaranteed not to be singular.
409
410 // Additionally, because Galois arithmetic is being used, all calulations
411 // involve exact values with no loss of precision. It is therefore
412 // not necessary to carry out any row or column swapping.
413
414 // Solve one row at a time
415
416 int progress = 0;
417
418 // For each row in the matrix
419 for (unsigned int row=0; row<datamissing; row++)
420 {
421 // NB Row and column swapping to find a non zero pivot value or to find the largest value
422 // is not necessary due to the nature of the arithmetic and construction of the RS matrix.
423
424 // Get the pivot value.
425 G pivotvalue = rightmatrix[row * rows + row];
426 assert(pivotvalue != 0);
427 if (pivotvalue == 0)
428 {
429 serr << "RS computation error." << endl;
430 return false;
431 }
432
433 // If the pivot value is not 1, then the whole row has to be scaled
434 if (pivotvalue != 1)
435 {
436 for (unsigned int col=0; col<leftcols; col++)
437 {
438 if (leftmatrix[row * leftcols + col] != 0)
439 {
440 leftmatrix[row * leftcols + col] /= pivotvalue;
441 }
442 }
443 rightmatrix[row * rows + row] = 1;
444 for (unsigned int col=row+1; col<rows; col++)
445 {
446 if (rightmatrix[row * rows + col] != 0)
447 {
448 rightmatrix[row * rows + col] /= pivotvalue;
449 }
450 }
451 }
452
453 // For every other row in the matrix
454 for (unsigned int row2=0; row2<rows; row2++)
455 {
456 // Define MPDL to skip reporting and speed things up
457 #ifndef MPDL
458 if (noiselevel > nlQuiet)
459 {
460 int newprogress = (row*rows+row2) * 1000 / (datamissing*rows);
461 if (progress != newprogress)
462 {
463 progress = newprogress;
464 sout << "Solving: " << progress/10 << '.' << progress%10 << "%\r" << flush;
465 }
466 }
467 #endif
468
469 if (row != row2)
470 {
471 // Get the scaling factor for this row.
472 G scalevalue = rightmatrix[row2 * rows + row];
473
474 if (scalevalue == 1)
475 {
476 // If the scaling factor happens to be 1, just subtract rows
477 for (unsigned int col=0; col<leftcols; col++)
478 {
479 if (leftmatrix[row * leftcols + col] != 0)
480 {
481 leftmatrix[row2 * leftcols + col] -= leftmatrix[row * leftcols + col];
482 }
483 }
484
485 for (unsigned int col=row; col<rows; col++)
486 {
487 if (rightmatrix[row * rows + col] != 0)
488 {
489 rightmatrix[row2 * rows + col] -= rightmatrix[row * rows + col];
490 }
491 }
492 }
493 else if (scalevalue != 0)
494 {
495 // If the scaling factor is not 0, then compute accordingly.
496 for (unsigned int col=0; col<leftcols; col++)
497 {
498 if (leftmatrix[row * leftcols + col] != 0)
499 {
500 leftmatrix[row2 * leftcols + col] -= leftmatrix[row * leftcols + col] * scalevalue;
501 }
502 }
503
504 for (unsigned int col=row; col<rows; col++)
505 {
506 if (rightmatrix[row * rows + col] != 0)
507 {
508 rightmatrix[row2 * rows + col] -= rightmatrix[row * rows + col] * scalevalue;
509 }
510 }
511 }
512 }
513 }
514 }
515 if (noiselevel > nlQuiet)
516 sout << "Solving: done." << endl;
517 if (noiselevel >= nlDebug)
518 {
519 for (unsigned int row=0; row<rows; row++)
520 {
521 sout << ((row==0) ? "/" : (row==rows-1) ? "\\" : "|");
522 for (unsigned int col=0; col<leftcols; col++)
523 {
524 sout << " "
525 << hex << setw(G::Bits>8?4:2) << setfill('0')
526 << (unsigned int)leftmatrix[row*leftcols+col];
527 }
528 sout << ((row==0) ? " \\ /" : (row==rows-1) ? " / \\" : " | |");
529 for (unsigned int col=0; col<rows; col++)
530 {
531 sout << " "
532 << hex << setw(G::Bits>8?4:2) << setfill('0')
533 << (unsigned int)rightmatrix[row*rows+col];
534 }
535 sout << ((row==0) ? " \\" : (row==rows-1) ? " /" : " | |");
536 sout << endl;
537
538 sout << dec << setw(0) << setfill(' ');
539 }
540 }
541
542 return true;
543 }
544
545
546 #endif // __REEDSOLOMON_H__
547