-
Notifications
You must be signed in to change notification settings - Fork 276
/
Copy pathRS.h
1107 lines (946 loc) · 50.6 KB
/
RS.h
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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*
* Copyright (C) 2014 Hard Consulting Corporation
* Copyright (C) 2023 Bryan Biedenkapp, N2PLL
* Copyright (C) 2024 Jonathan Naylor, G4KLX
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
#if !defined(RS_H)
#define RS_H
/*
* Ezpwd Reed-Solomon -- Reed-Solomon encoder / decoder library
*
* The core Reed-Solomon codec implementation in c++/ezpwd/rs_base is by Phil Karn, converted to C++
* by Perry Kundert ([email protected]), and may be used under the terms of the LGPL. Here
* is the terms from Phil's README file (see phil-karn/fec-3.0.1/README):
*
* COPYRIGHT
*
* This package is copyright 2006 by Phil Karn, KA9Q. It may be used
* under the terms of the GNU Lesser General Public License (LGPL). See
* the file "lesser.txt" in this package for license details.
*
* The c++/ezpwd/rs_base file is, therefore, redistributed under the terms of the LGPL, while the
* rest of Ezpwd Reed-Solomon is distributed under either the GPL or Commercial licenses.
* Therefore, even if you have obtained Ezpwd Reed-Solomon under a Commercial license, you must make
* available the source code of the c++/ezpwd/rs_base file with your product. One simple way to
* accomplish this is to include the following URL in your code or documentation:
*
* https://github.com/pjkundert/ezpwd-reed-solomon/blob/master/c++/ezpwd/rs_base
*
*
* The Linux 3.15.1 version of lib/reed_solomon was also consulted as a cross-reference, which (in
* turn) is basically verbatim copied from Phil Karn's LGPL implementation, to ensure that no new
* defects had been found and fixed; there were no meaningful changes made to Phil's implementation.
* I've personally been using Phil's implementation for years in a heavy industrial use, and it is
* rock-solid.
*
* However, both Phil's and the Linux kernel's (copy of Phil's) implementation will return a
* "corrected" decoding with impossible error positions, in some cases where the error load
* completely overwhelms the R-S encoding. These cases, when detected, are rejected in this
* implementation. This could be considered a defect in Phil's (and hence the Linux kernel's)
* implementations, which results in them accepting clearly incorrect R-S decoded values as valid
* (corrected) R-S codewords. We chose the report failure on these attempts.
*
*/
#include <algorithm>
#include <array>
#include <cstdint>
#include <cstring>
#include <iostream>
#include <iomanip>
#include <type_traits>
#include <vector>
//
// Preprocessor defines available:
//
// EZPWD_NO_EXCEPTS -- define to use no exceptions; return -1, or abort on catastrophic failures
// EZPWD_NO_MOD_TAB -- define to force no "modnn" Galois modulo table acceleration
// EZPWD_ARRAY_SAFE -- define to force usage of bounds-checked arrays for most tabular data
// EZPWD_ARRAY_TEST -- define to force erroneous sizing of some arrays for non-production testing
//
#if defined(EZPWD_NO_EXCEPTS)
#include <cstdio>
#define EZPWD_RAISE_OR_ABORT(typ, str) do { \
std::fputs((str), stderr); std::fputc('\n', stderr);\
abort(); \
} while (false)
#define EZPWD_RAISE_OR_RETURN(typ, str, ret) return (ret)
#else
#define EZPWD_RAISE_OR_ABORT(typ, str) throw (typ)(str)
#define EZPWD_RAISE_OR_RETURN(typ, str, ret) throw (typ)(str)
#endif
namespace rs
{
// ezpwd::log_<N,B> -- compute the log base B of N at compile-time
template <size_t N, size_t B = 2> struct log_ { enum { value = 1 + log_<N / B, B>::value }; };
template <size_t B> struct log_<1, B> { enum { value = 0 }; };
template <size_t B> struct log_<0, B> { enum { value = 0 }; };
// ---------------------------------------------------------------------------
// Class Declaration
// ---------------------------------------------------------------------------
/**
* @brief Reed-Solomon codec generic base class.
* @ingroup edac_rs
*/
class reed_solomon_base {
public:
/** @brief A data element's bits. */
virtual size_t datum() const = 0;
/** @brief A symbol's bits. */
virtual size_t symbol() const = 0;
/** @brief R-S block size (maximum total symbols). */
virtual int size() const = 0;
/** @brief R-S roots (parity symbols). */
virtual int nroots() const = 0;
/** @brief R-S net payload (data symbols). */
virtual int load() const = 0;
/** @brief Initializes a new instance of the reed_solomon_base class. */
reed_solomon_base() = default;
/** @brief Finalizes a instance of the reed_solomon_base class. */
virtual ~reed_solomon_base() = default;
/** @brief */
virtual std::ostream& output(std::ostream& lhs) const { return lhs << "RS(" << this->size() << "," << this->load() << ")"; }
//
// {en,de}code -- Compute/Correct errors/erasures in a Reed-Solomon encoded container
//
// The encoded parity symbols may be included in 'data' (len includes nroots() parity
// symbols), or may (optionally) supplied separately in (at least nroots()-sized)
// 'parity'.
//
// For decode, optionally specify some known erasure positions (up to nroots()). If
// non-empty 'erasures' is provided, it contains the positions of each erasure. If a
// non-zero pointer to a 'position' vector is provided, its capacity will be increased to
// be capable of storing up to 'nroots()' ints; the actual deduced error locations will be
// returned.
//
// RETURN VALUE
//
// Return -1 on error. The encode returns the number of parity symbols produced;
// decode returns the number of symbols corrected. Both errors and erasures are included,
// so long as they are actually different than the deduced value. In other words, if a
// symbol is marked as an erasure but it actually turns out to be correct, it's index will
// NOT be included in the returned count, nor the modified erasure vector!
//
int encode(std::string& data) const
{
typedef uint8_t uT;
typedef std::pair<uT*, uT*> uTpair;
data.resize(data.size() + nroots());
return encode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()));
}
int encode(const std::string& data, std::string& parity) const
{
typedef uint8_t uT;
typedef std::pair<const uT*, const uT*> cuTpair;
typedef std::pair<uT*, uT*> uTpair;
parity.resize(nroots());
return encode(cuTpair((const uT*)&data.front(), (const uT*)&data.front() + data.size()),
uTpair((uT*)&parity.front(), (uT*)&parity.front() + parity.size()));
}
template<typename T> int encode(std::vector<T>& data) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<uT*, uT*> uTpair;
data.resize(data.size() + nroots());
return encode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()));
}
template<typename T> int encode(const std::vector<T>& data, std::vector<T>& parity) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<const uT*, const uT*> cuTpair;
typedef std::pair<uT*, uT*> uTpair;
parity.resize(nroots());
return encode(cuTpair((uT*)&data.front(), (uT*)&data.front() + data.size()),
uTpair((uT*)&parity.front(), (uT*)&parity.front() + parity.size()));
}
template<typename T, size_t N> int encode(std::array<T, N>& data, int pad = 0) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<uT*, uT*> uTpair;
return encode(uTpair((uT*)&data.front() + pad, (uT*)&data.front() + data.size()));
}
virtual int encode(const std::pair<uint8_t*, uint8_t*>& data) const = 0;
virtual int encode(const std::pair<const uint8_t*, const uint8_t*>& data,
const std::pair<uint8_t*, uint8_t*>& parity) const = 0;
virtual int encode(const std::pair<uint16_t*, uint16_t*>& data) const = 0;
virtual int encode(const std::pair<const uint16_t*, const uint16_t*>& data,
const std::pair<uint16_t*, uint16_t*>& parity) const = 0;
virtual int encode(const std::pair<uint32_t*, uint32_t*>& data) const = 0;
virtual int encode(const std::pair<const uint32_t*, const uint32_t*>& data,
const std::pair<uint32_t*, uint32_t*>& parity) const = 0;
int decode(std::string& data, const std::vector<int>& erasure = std::vector<int>(),
std::vector<int>* position = 0) const
{
typedef uint8_t uT;
typedef std::pair<uT*, uT*> uTpair;
return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()), erasure, position);
}
int decode(std::string& data, std::string& parity, const std::vector<int>& erasure = std::vector<int>(),
std::vector<int>* position = 0) const
{
typedef uint8_t uT;
typedef std::pair<uT*, uT*> uTpair;
return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()),
uTpair((uT*)&parity.front(), (uT*)&parity.front() + parity.size()), erasure, position);
}
template<typename T> int decode(std::vector<T>& data, const std::vector<int>& erasure = std::vector<int>(),
std::vector<int>* position = 0) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<uT*, uT*> uTpair;
return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()), erasure, position);
}
template<typename T> int decode(std::vector<T>& data, std::vector<T>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<uT*, uT*> uTpair;
return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()),
uTpair((uT*)&parity.front(), (uT*)&parity.front() + parity.size()), erasure, position);
}
template<typename T, size_t N> int decode(std::array<T, N>& data, int pad = 0,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<uT*, uT*> uTpair;
return decode(uTpair((uT*)&data.front() + pad, (uT*)&data.front() + data.size()), erasure, position);
}
virtual int decode(const std::pair<uint8_t*, uint8_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
virtual int decode(const std::pair<uint8_t*, uint8_t*>& data, const std::pair<uint8_t*, uint8_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
virtual int decode(const std::pair<uint16_t*, uint16_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
virtual int decode(const std::pair<uint16_t*, uint16_t*>& data, const std::pair<uint16_t*, uint16_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
virtual int decode(const std::pair<uint32_t*, uint32_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
virtual int decode(const std::pair<uint32_t*, uint32_t*>& data, const std::pair<uint32_t*, uint32_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
};
//
// std::ostream << edac::rs::reed_solomon<...>
//
// Output a R-S codec description in standard form eg. RS(255,253)
//
inline std::ostream& operator<<(std::ostream& lhs, const rs::reed_solomon_base& rhs) { return rhs.output(lhs); }
// ---------------------------------------------------------------------------
// Structure Declaration
// ---------------------------------------------------------------------------
/**
* @brief Default field polynomial generator functor.
* @tparam SYM
* @tparam PLY
* @ingroup edac_rs
*/
template<int SYM, int PLY>
struct gfpoly {
int operator() (int sr) const
{
if (sr == 0) {
sr = 1;
} else {
sr <<= 1;
if (sr & (1 << SYM))
sr ^= PLY;
sr &= ((1 << SYM) - 1);
}
return sr;
}
};
// ---------------------------------------------------------------------------
// Class Declaration
// ---------------------------------------------------------------------------
/**
* @brief R-S tables common to all RS(NN,*) with same SYM, PRM and PLY.
* @tparam TYP
* @tparam SYM
* @tparam PRM
* @tparam PLY
* @ingroup edac_rs
*/
template <typename TYP, int SYM, int PRM, class PLY>
class reed_solomon_tabs : public reed_solomon_base {
public:
typedef TYP symbol_t;
/** @brief Bits / TYP */
static const size_t DATUM = 8 * sizeof TYP();
/** @brief Bits / Symbol */
static const size_t SYMBOL = SYM;
static const int MM = SYM;
static const int SIZE = (1 << SYM) - 1; // maximum symbols in field
static const int NN = SIZE;
static const int A0 = SIZE;
// modulo table: 1/2 the symbol size squared, up to 4k
static const int MODS = SYM > 8 ? (1 << 12) : (1 << SYM << SYM / 2);
static int iprim; // initialized to -1, below
protected:
static std::array<TYP, NN + 1> alpha_to;
static std::array<TYP, NN + 1> index_of;
static std::array<TYP, MODS> mod_of;
/** @brief Initializes a new instance of the reed_solomon_tabs class. */
reed_solomon_tabs() : reed_solomon_base()
{
// Do init if not already done. We check one value which is initialized to -1; this is
// safe, 'cause the value will not be set 'til the initializing thread has completely
// initialized the structure. Worst case scenario: multiple threads will initialize
// identically. No mutex necessary.
if (iprim >= 0)
return;
#if DEBUG_RS
LogDebug(LOG_HOST, "reed_solomon_tabs::reed_solomon_tabs() RS(%d,*): initialized for %d symbols size, %d modulo table", SIZE, NN, MODS);
#endif
// Generate Galois field lookup tables
index_of[0] = A0; // log(zero) = -inf
alpha_to[A0] = 0; // alpha**-inf = 0
PLY poly;
int sr = poly(0);
for (int i = 0; i < NN; i++) {
index_of[sr] = i;
alpha_to[i] = sr;
sr = poly(sr);
}
// If it's not primitive, raise exception or abort
if (sr != alpha_to[0])
EZPWD_RAISE_OR_ABORT(std::runtime_error, "reed-solomon: Galois field polynomial not primitive");
// Generate modulo table for some commonly used (non-trivial) values
for (int x = NN; x < NN + MODS; ++x)
mod_of[x - NN] = _modnn(x);
// Find prim-th root of 1, index form, used in decoding.
int iptmp = 1;
while (iptmp % PRM != 0)
iptmp += NN;
iprim = iptmp / PRM;
}
/// <summary>Finalizes a instance of the reed_solomon_tabs class.</summary>
~reed_solomon_tabs() override = default;
//
// modnn -- modulo replacement for galois field arithmetics, optionally w/ table acceleration
//
// @x: the value to reduce (will never be -'ve)
//
// where
// MM = number of bits per symbol
// NN = (2^MM) - 1
//
// Simple arithmetic modulo would return a wrong result for values >= 3 * NN
//
TYP _modnn(int x) const
{
while (x >= NN) {
x -= NN;
x = (x >> MM) + (x & NN);
}
return x;
}
TYP modnn(int x) const
{
while (x >= NN + MODS) {
x -= NN;
x = (x >> MM) + (x & NN);
}
if (MODS && x >= NN)
x = mod_of[x - NN];
return x;
}
};
// ---------------------------------------------------------------------------
// Class Declaration
// ---------------------------------------------------------------------------
/**
* @brief Reed-Solomon codec.
* @tparam TYP A symbol datum; {en,de}code operates on arrays of these
* @tparam SYM Bits per symbol
* @tparam RTS
* @tparam FCR First consecutive root, index form
* @tparam PRM Primitive element, index form
* @tparam PLY The primitive generator polynominal functor
* @ingroup edac_rs
*/
/*
* @TYP: A symbol datum; {en,de}code operates on arrays of these
* @DATUM: Bits per datum (a TYP())
* @SYM{BOL}, MM: Bits per symbol
* @NN: Symbols per block (== (1<<MM)-1)
* @alpha_to: log lookup table
* @index_of: Antilog lookup table
* @genpoly: Generator polynomial
* @NROOTS: Number of generator roots = number of parity symbols
* @FCR: First consecutive root, index form
* @PRM: Primitive element, index form
* @iprim: prim-th root of 1, index form
* @PLY: The primitive generator polynominal functor
*
* All reed_solomon<T, ...> instances with the same template type parameters share a common
* (static) set of alpha_to, index_of and genpoly tables. The first instance to be constructed
* initializes the tables.
*
* Each specialized type of reed_solomon implements a specific encode/decode method
* appropriate to its datum 'TYP'. When accessed via a generic reed_solomon_base pointer, only
* access via "safe" (size specifying) containers or iterators is available.
*/
template<typename TYP, int SYM, int RTS, int FCR, int PRM, class PLY>
class reed_solomon : public reed_solomon_tabs<TYP, SYM, PRM, PLY> {
public:
typedef reed_solomon_tabs<TYP, SYM, PRM, PLY> tabs_t;
using tabs_t::DATUM;
using tabs_t::SYMBOL;
using tabs_t::MM;
using tabs_t::SIZE;
using tabs_t::NN;
using tabs_t::A0;
using tabs_t::iprim;
using tabs_t::alpha_to;
using tabs_t::index_of;
using tabs_t::modnn;
static const int NROOTS = RTS;
static const int LOAD = SIZE - NROOTS; // maximum non-parity symbol payload
protected:
static std::array<TYP, NROOTS + 1> genpoly;
public:
/** @brief Initializes a new instance of the reed_solomon class. */
reed_solomon() : reed_solomon_tabs<TYP, SYM, PRM, PLY>()
{
// We check one element of the array; this is safe, 'cause the value will not be
// initialized 'til the initializing thread has completely initialized the array. Worst
// case scenario: multiple threads will initialize identically. No mutex necessary.
if (genpoly[0])
return;
#if DEBUG_RS
LogDebug(LOG_HOST, "reed_solomon::reed_solomon() RS(%d,%d): initialized for %d roots", SIZE, LOAD, NROOTS);
#endif
std::array<TYP, NROOTS + 1> tmppoly; // uninitialized
// Form RS code generator polynomial from its roots. Only lower-index entries are
// consulted, when computing subsequent entries; only index 0 needs initialization.
tmppoly[0] = 1;
for (int i = 0, root = FCR * PRM; i < NROOTS; i++, root += PRM) {
tmppoly[i + 1] = 1;
// Multiply tmppoly[] by @**(root + x)
for (int j = i; j > 0; j--) {
if (tmppoly[j] != 0)
tmppoly[j] = tmppoly[j - 1] ^ alpha_to[modnn(index_of[tmppoly[j]] + root)];
else
tmppoly[j] = tmppoly[j - 1];
}
// tmppoly[0] can never be zero
tmppoly[0] = alpha_to[modnn(index_of[tmppoly[0]] + root)];
}
// convert NROOTS entries of tmppoly[] to genpoly[] in index form for quicker encoding,
// in reverse order so genpoly[0] is last element initialized.
for (int i = NROOTS; i >= 0; --i)
genpoly[i] = index_of[tmppoly[i]];
}
/** @brief Finalizes a instance of the reed_solomon class. */
virtual ~reed_solomon() = default;
/** @brief A data element's bits. */
virtual size_t datum() const { return DATUM; }
/** @brief A symbol's bits. */
virtual size_t symbol() const { return SYMBOL; }
/** @brief R-S block size (maximum total symbols). */
virtual int size() const { return SIZE; }
/** @brief R-S roots (parity symbols). */
virtual int nroots() const { return NROOTS; }
/** @brief R-S net payload (data symbols). */
virtual int load() const { return LOAD; }
using reed_solomon_base::encode;
virtual int encode(const std::pair<uint8_t*, uint8_t*>& data) const { return encode_mask(data.first, int(data.second - data.first) - NROOTS, data.second - NROOTS); }
virtual int encode(const std::pair<const uint8_t*, const uint8_t*>& data,
const std::pair<uint8_t*, uint8_t*>& parity) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return encode_mask(data.first, int(data.second - data.first), parity.first);
}
virtual int encode(const std::pair<uint16_t*, uint16_t*>& data) const { return encode_mask(data.first, int(data.second - data.first) - NROOTS, data.second - NROOTS); }
virtual int encode(const std::pair<const uint16_t*, const uint16_t*>& data,
const std::pair<uint16_t*, uint16_t*>& parity) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return encode_mask(data.first, int(data.second - data.first), parity.first);
}
virtual int encode(const std::pair<uint32_t*, uint32_t*>& data) const { return encode_mask(data.first, int(data.second - data.first) - NROOTS, data.second - NROOTS); }
virtual int encode(const std::pair<const uint32_t*, const uint32_t*>& data,
const std::pair<uint32_t*, uint32_t*>& parity) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return encode_mask(data.first, int(data.second - data.first), parity.first);
}
template<typename INP>
int encode_mask(const INP* data, int len, INP* parity) const
{
if (len < 1)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: must provide space for all parity and at least one non-parity symbol", -1);
const TYP* dataptr;
TYP* pariptr;
const size_t INPUT = 8 * sizeof(INP);
if (DATUM != SYMBOL || DATUM != INPUT) {
// Our DATUM (TYP) size (eg. uint8_t ==> 8, uint16_t ==> 16, uint32_t ==> 32)
// doesn't exactly match our R-S SYMBOL size (eg. 6), or our INP size; Must mask and
// copy. The INP data must fit at least the SYMBOL size!
if (SYMBOL > INPUT)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: output data type too small to contain symbols", -1);
std::array<TYP, SIZE> tmp;
TYP msk = static_cast<TYP>(~0UL << SYMBOL);
for (int i = 0; i < len; ++i)
tmp[LOAD - len + i] = data[i] & ~msk;
dataptr = &tmp[LOAD - len];
pariptr = &tmp[LOAD];
encode(dataptr, len, pariptr);
// we copied/masked data; copy the parity symbols back (may be different sizes)
for (int i = 0; i < NROOTS; ++i)
parity[i] = pariptr[i];
} else {
// Our R-S SYMBOL size, DATUM size and INP type size exactly matches; use in-place.
dataptr = reinterpret_cast<const TYP*>(data);
pariptr = reinterpret_cast<TYP*>(parity);
encode(dataptr, len, pariptr);
}
return NROOTS;
}
using reed_solomon_base::decode;
virtual int decode(const std::pair<uint8_t*, uint8_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
return decode_mask(data.first, int(data.second - data.first), (uint8_t*)0, erasure, position);
}
virtual int decode(const std::pair<uint8_t*, uint8_t*>& data, const std::pair<uint8_t*, uint8_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return decode_mask(data.first, int(data.second - data.first), parity.first, erasure, position);
}
virtual int decode(const std::pair<uint16_t*, uint16_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
return decode_mask(data.first, int(data.second - data.first), (uint16_t*)0, erasure, position);
}
virtual int decode(const std::pair<uint16_t*, uint16_t*>& data, const std::pair<uint16_t*, uint16_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return decode_mask(data.first, int(data.second - data.first), parity.first, erasure, position);
}
virtual int decode(const std::pair<uint32_t*, uint32_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
return decode_mask(data.first, int(data.second - data.first), (uint32_t*)0, erasure, position);
}
virtual int decode(const std::pair<uint32_t*, uint32_t*>& data, const std::pair<uint32_t*, uint32_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return decode_mask(data.first, int(data.second - data.first), parity.first, erasure, position);
}
//
// decode_mask -- mask INP data into valid SYMBOL data
//
// Incoming data may be in a variety of sizes, and may contain information beyond the
// R-S symbol capacity. For example, we might use a 6-bit R-S symbol to correct the lower
// 6 bits of an 8-bit data character. This would allow us to correct common substitution
// errors (such as '2' for '3', 'R' for 'T', 'n' for 'm').
//
template<typename INP>
int decode_mask(INP* data, int len, INP* parity = 0, const std::vector<int>& erasure = std::vector<int>(),
std::vector<int>* position = 0) const {
if (len < (parity ? 0 : NROOTS) + 1)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: must provide all parity and at least one non-parity symbol", -1);
if (!parity) {
len -= NROOTS;
parity = data + len;
}
TYP* dataptr;
TYP* pariptr;
const size_t INPUT = 8 * sizeof(INP);
std::array<TYP, SIZE> tmp;
TYP msk = static_cast<TYP>(~0UL << SYMBOL);
const bool cpy = DATUM != SYMBOL || DATUM != INPUT;
if (cpy) {
// Our DATUM (TYP) size (eg. uint8_t ==> 8, uint16_t ==> 16, uint32_t ==> 32)
// doesn't exactly match our R-S SYMBOL size (eg. 6), or our INP size; Must copy.
// The INP data must fit at least the SYMBOL size!
if (SYMBOL > INPUT)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: input data type too small to contain symbols", -1);
for (int i = 0; i < len; ++i)
tmp[LOAD - len + i] = data[i] & ~msk;
dataptr = &tmp[LOAD - len];
for (int i = 0; i < NROOTS; ++i) {
if (TYP(parity[i]) & msk)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity data contains information beyond R-S symbol size", -1);
tmp[LOAD + i] = (TYP)parity[i];
}
pariptr = &tmp[LOAD];
} else {
// Our R-S SYMBOL size, DATUM size and INPUT type sizes exactly matches
dataptr = reinterpret_cast<TYP*>(data);
pariptr = reinterpret_cast<TYP*>(parity);
}
int corrects;
if (!erasure.size() && !position) {
// No erasures, and error position info not wanted.
corrects = decode(dataptr, len, pariptr);
} else {
// Either erasure location info specified, or resultant error position info wanted;
// Prepare pos (a temporary, if no position vector provided), and copy any provided
// erasure positions. After number of corrections is known, resize the position
// vector. Thus, we use any supplied erasure info, and optionally return any
// correction position info separately.
std::vector<int> _pos;
std::vector<int>& pos = position ? *position : _pos;
pos.resize(std::max(size_t(NROOTS), erasure.size()));
std::copy(erasure.begin(), erasure.end(), pos.begin());
corrects = decode(dataptr, len, pariptr, &pos.front(), int(erasure.size()));
if (corrects > int(pos.size()))
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: FATAL: produced too many corrections; possible corruption!", -1);
pos.resize(std::max(0, corrects));
}
if (cpy && corrects > 0) {
for (int i = 0; i < len; ++i) {
data[i] &= msk;
data[i] |= tmp[LOAD - len + i];
}
for (int i = 0; i < NROOTS; ++i)
parity[i] = tmp[LOAD + i];
}
return corrects;
}
int encode(const TYP* data, int len, TYP* parity) const
{
// Check length parameter for validity
int pad = NN - NROOTS - len;
if (pad < 0 || pad >= NN)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: data length incompatible with block size and error correction symbols", -1);
for (int i = 0; i < NROOTS; i++)
parity[i] = 0;
for (int i = 0; i < len; i++) {
TYP feedback = index_of[data[i] ^ parity[0]];
if (feedback != A0) {
for (int j = 1; j < NROOTS; j++)
parity[j] ^= alpha_to[modnn(feedback + genpoly[NROOTS - j])];
}
std::rotate(parity, parity + 1, parity + NROOTS);
if (feedback != A0)
parity[NROOTS - 1] = alpha_to[modnn(feedback + genpoly[0])];
else
parity[NROOTS - 1] = 0;
}
return NROOTS;
}
int decode(TYP* data, int len, TYP* parity, int* eras_pos = 0, int no_eras = 0, TYP* corr = 0) const
{
typedef std::array<TYP, NROOTS> typ_nroots;
typedef std::array<TYP, NROOTS + 1> typ_nroots_1;
typedef std::array<int, NROOTS> int_nroots;
typ_nroots_1 lambda{ {0} };
typ_nroots syn;
typ_nroots_1 b;
typ_nroots_1 t;
typ_nroots_1 omega;
int_nroots root;
typ_nroots_1 reg;
int_nroots loc;
int count = 0;
// Check length parameter and erasures for validity
int pad = NN - NROOTS - len;
if (pad < 0 || pad >= NN)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: data length incompatible with block size and error correction symbols", -1);
if (no_eras) {
if (no_eras > NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: number of erasures exceeds capacity (number of roots)", -1);
for (int i = 0; i < no_eras; ++i) {
if (eras_pos[i] < 0 || eras_pos[i] >= len + NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: erasure positions outside data+parity", -1);
}
}
// form the syndromes; i.e., evaluate data(x) at roots of g(x)
for (int i = 0; i < NROOTS; i++)
syn[i] = data[0];
for (int j = 1; j < len; j++) {
for (int i = 0; i < NROOTS; i++) {
if (syn[i] == 0)
syn[i] = data[j];
else
syn[i] = data[j] ^ alpha_to[modnn(index_of[syn[i]] + (FCR + i) * PRM)];
}
}
for (int j = 0; j < NROOTS; j++) {
for (int i = 0; i < NROOTS; i++) {
if (syn[i] == 0)
syn[i] = parity[j];
else
syn[i] = parity[j] ^ alpha_to[modnn(index_of[syn[i]] + (FCR + i) * PRM)];
}
}
// Convert syndromes to index form, checking for nonzero condition
TYP syn_error = 0;
for (int i = 0; i < NROOTS; i++) {
syn_error |= syn[i];
syn[i] = index_of[syn[i]];
}
int deg_lambda = 0;
int deg_omega = 0;
int r = no_eras;
int el = no_eras;
if (!syn_error) {
// if syndrome is zero, data[] is a codeword and there are no errors to correct.
count = 0;
goto finish; // ewww; gotos!
}
lambda[0] = 1;
if (no_eras > 0) {
// Init lambda to be the erasure locator polynomial. Convert erasure positions
// from index into data, to index into Reed-Solomon block.
lambda[1] = alpha_to[modnn(PRM * (NN - 1 - (eras_pos[0] + pad)))];
for (int i = 1; i < no_eras; i++) {
TYP u = modnn(PRM * (NN - 1 - (eras_pos[i] + pad)));
for (int j = i + 1; j > 0; j--) {
TYP tmp = index_of[lambda[j - 1]];
if (tmp != A0)
lambda[j] ^= alpha_to[modnn(u + tmp)];
}
}
}
#if DEBUG_RS
// Test code that verifies the erasure locator polynomial just constructed
// Needed only for decoder debugging.
// find roots of the erasure location polynomial
for (int i = 1; i <= no_eras; i++) {
reg[i] = index_of[lambda[i]];
}
count = 0;
for (int i = 1, k = iprim - 1; i <= NN; i++, k = modnn(k + iprim)) {
TYP q = 1;
for (int j = 1; j <= no_eras; j++) {
if (reg[j] != A0) {
reg[j] = modnn(reg[j] + j);
q ^= alpha_to[reg[j]];
}
}
if (q != 0) {
continue;
}
// store root and error location number indices
root[count] = i;
loc[count] = k;
count++;
}
if (count != no_eras) {
LogDebug(LOG_HOST, "reed_solomon::decode(): count = %d, no_eras = %d, lambda(x) is WRONG", count, no_eras);
count = -1;
goto finish;
}
if (count) {
std::stringstream ss;
ss << "reed_solomon::decode(): Erasure positions as determined by roots of Eras Loc Poly: ";
for (int i = 0; i < count; i++) {
ss << loc[i] << ' ';
}
LogDebug(LOG_HOST, "%s", ss.str().c_str());
ss.clear();
ss << "reed_solomon::decode(): Erasure positions as determined by roots of eras_pos array: ";
for (int i = 0; i < no_eras; i++) {
ss << eras_pos[i] << ' ';
}
LogDebug(LOG_HOST, "%s", ss.str().c_str());
}
#endif
for (int i = 0; i < NROOTS + 1; i++)
b[i] = index_of[lambda[i]];
//
// Begin Berlekamp-Massey algorithm to determine error+erasure locator polynomial
//
while (++r <= NROOTS) {
// r is the step number
// Compute discrepancy at the r-th step in poly-form
TYP discr_r = 0;
for (int i = 0; i < r; i++) {
if ((lambda[i] != 0) && (syn[r - i - 1] != A0))
discr_r ^= alpha_to[modnn(index_of[lambda[i]] + syn[r - i - 1])];
}
discr_r = index_of[discr_r]; // Index form
if (discr_r == A0) {
// 2 lines below: B(x) <-- x*B(x)
// Rotate the last element of b[NROOTS+1] to b[0]
std::rotate(b.begin(), b.begin() + NROOTS, b.end());
b[0] = A0;
} else {
// 7 lines below: T(x) <-- lambda(x)-discr_r*x*b(x)
t[0] = lambda[0];
for (int i = 0; i < NROOTS; i++) {
if (b[i] != A0)
t[i + 1] = lambda[i + 1] ^ alpha_to[modnn(discr_r + b[i])];
else
t[i + 1] = lambda[i + 1];
}
if (2 * el <= r + no_eras - 1) {
el = r + no_eras - el;
// 2 lines below: B(x) <-- inv(discr_r) * lambda(x)
for (int i = 0; i <= NROOTS; i++)
b[i] = ((lambda[i] == 0) ? A0 : modnn(index_of[lambda[i]] - discr_r + NN));
} else {
// 2 lines below: B(x) <-- x*B(x)
std::rotate(b.begin(), b.begin() + NROOTS, b.end());
b[0] = A0;
}
lambda = t;
}
}
// Convert lambda to index form and compute deg(lambda(x))
for (int i = 0; i < NROOTS + 1; i++) {
lambda[i] = index_of[lambda[i]];
if (lambda[i] != NN)
deg_lambda = i;
}
// Find roots of error+erasure locator polynomial by Chien search
reg = lambda;
count = 0; // Number of roots of lambda(x)
for (int i = 1, k = iprim - 1; i <= NN; i++, k = modnn(k + iprim)) {
TYP q = 1; // lambda[0] is always 0
for (int j = deg_lambda; j > 0; j--) {
if (reg[j] != A0) {
reg[j] = modnn(reg[j] + j);
q ^= alpha_to[reg[j]];
}
}
if (q != 0)
continue; // Not a root
// store root (index-form) and error location number
#if DEBUG_RS
LogDebug(LOG_HOST, "reed_solomon::decode(): count = %d, root = %d, loc = %d", count, i, k);
#endif
root[count] = i;
loc[count] = k;
// If we've already found max possible roots, abort the search to save time
if (++count == deg_lambda)
break;
}
if (deg_lambda != count) {
// deg(lambda) unequal to number of roots => uncorrectable error detected
count = -1;
goto finish;
}
//
// Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo x**NROOTS). in
// index form. Also find deg(omega).
//
deg_omega = deg_lambda - 1;
for (int i = 0; i <= deg_omega; i++) {
TYP tmp = 0;
for (int j = i; j >= 0; j--) {
if ((syn[i - j] != A0) && (lambda[j] != A0)) {
tmp ^= alpha_to[modnn(syn[i - j] + lambda[j])];
}
}
omega[i] = index_of[tmp];
}
//
// Compute error values in poly-form. num1 = omega(inv(X(l))), num2 = inv(X(l))**(fcr-1)