用于 Reed-Solomon 解码的勘误表（擦除+错误）Berlekamp-Massey
Errata (erasures+errors) Berlekamp-Massey for Reed-Solomon decoding
我正在尝试在 Python 中实现支持擦除和错误解码的 Reed-Solomon 编码器-解码器，这让我抓狂。
该实现目前支持仅解码错误或仅解码擦除，但不能同时解码（即使它低于 2*errors+erasures <= (n-k) 的理论界限）。
从 Blahut 的论文 (here and here) 看来，我们只需要用擦除定位器多项式初始化错误定位器多项式就可以隐式计算 Berlekamp-Massey 中的勘误定位器多项式。
这种方法对我部分有效：当我有 2*errors+erasures < (n-k)/2 时它有效，但实际上调试后它只有效，因为 BM 计算了一个错误定位器多项式，它得到了完全相同的值作为擦除定位器多项式（因为我们低于仅错误校正的限制），因此它通过伽罗瓦域被截断，我们最终得到擦除定位器多项式的正确值（至少我是这样理解的，我可能是错误的）。
然而，当我们超过 (n-k)/2 时，例如，如果 n = 20 和 k = 11，那么我们有 (n-k)=9 个我们可以纠正的擦除符号，如果我们输入 5 个擦除，那么 BM只是出错了。如果我们输入 4 个擦除 + 1 个错误（我们仍然远低于界限，因为我们有 2*错误+擦除 = 2+4 = 6 < 9），BM 仍然出错。
我实现的Berlekamp-Massey的具体算法可以在this presentation (pages 15-17), but a very similar description can be found here and here中找到，这里我附上一份数学描述：
现在，我几乎可以将这个数学算法精确地复制到 Python 代码中。我想要的是扩展它以支持擦除，我尝试通过使用擦除定位器初始化错误定位器 sigma:
def _berlekamp_massey(self, s, k=None, erasures_loc=None):
'''Computes and returns the error locator polynomial (sigma) and the
error evaluator polynomial (omega).
If the erasures locator is specified, we will return an errors-and-erasures locator polynomial and an errors-and-erasures evaluator polynomial.
The parameter s is the syndrome polynomial (syndromes encoded in a
generator function) as returned by _syndromes. Don't be confused with
the other s = (n-k)/2
Notes:
The error polynomial:
E(x) = E_0 + E_1 x + ... + E_(n-1) x^(n-1)
j_1, j_2, ..., j_s are the error positions. (There are at most s
errors)
Error location X_i is defined: X_i = a^(j_i)
that is, the power of a corresponding to the error location
Error magnitude Y_i is defined: E_(j_i)
that is, the coefficient in the error polynomial at position j_i
Error locator polynomial:
sigma(z) = Product( 1 - X_i * z, i=1..s )
roots are the reciprocals of the error locations
( 1/X_1, 1/X_2, ...)
Error evaluator polynomial omega(z) is here computed at the same time as sigma, but it can also be constructed afterwards using the syndrome and sigma (see _find_error_evaluator() method).
'''
# For errors-and-erasures decoding, see: Blahut, Richard E. "Transform techniques for error control codes." IBM Journal of Research and development 23.3 (1979): 299-315. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.600&rep=rep1&type=pdf and also a MatLab implementation here: http://www.mathworks.com/matlabcentral/fileexchange/23567-reed-solomon-errors-and-erasures-decoder/content//RS_E_E_DEC.m
# also see: Blahut, Richard E. "A universal Reed-Solomon decoder." IBM Journal of Research and Development 28.2 (1984): 150-158. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.84.2084&rep=rep1&type=pdf
# or alternatively see the reference book by Blahut: Blahut, Richard E. Theory and practice of error control codes. Addison-Wesley, 1983.
# and another good alternative book with concrete programming examples: Jiang, Yuan. A practical guide to error-control coding using Matlab. Artech House, 2010.
n = self.n
if not k: k = self.k
# Initialize:
if erasures_loc:
sigma = [ Polynomial(erasures_loc.coefficients) ] # copy erasures_loc by creating a new Polynomial
B = [ Polynomial(erasures_loc.coefficients) ]
else:
sigma = [ Polynomial([GF256int(1)]) ] # error locator polynomial. Also called Lambda in other notations.
B = [ Polynomial([GF256int(1)]) ] # this is the error locator support/secondary polynomial, which is a funky way to say that it's just a temporary variable that will help us construct sigma, the error locator polynomial
omega = [ Polynomial([GF256int(1)]) ] # error evaluator polynomial. We don't need to initialize it with erasures_loc, it will still work, because Delta is computed using sigma, which itself is correctly initialized with erasures if needed.
A = [ Polynomial([GF256int(0)]) ] # this is the error evaluator support/secondary polynomial, to help us construct omega
L = [ 0 ] # necessary variable to check bounds (to avoid wrongly eliminating the higher order terms). For more infos, see https://www.cs.duke.edu/courses/spring11/cps296.3/decoding_rs.pdf
M = [ 0 ] # optional variable to check bounds (so that we do not mistakenly overwrite the higher order terms). This is not necessary, it's only an additional safe check. For more infos, see the presentation decoding_rs.pdf by Andrew Brown in the doc folder.
# Polynomial constants:
ONE = Polynomial(z0=GF256int(1))
ZERO = Polynomial(z0=GF256int(0))
Z = Polynomial(z1=GF256int(1)) # used to shift polynomials, simply multiply your poly * Z to shift
s2 = ONE + s
# Iteratively compute the polynomials 2s times. The last ones will be
# correct
for l in xrange(0, n-k):
K = l+1
# Goal for each iteration: Compute sigma[K] and omega[K] such that
# (1 + s)*sigma[l] == omega[l] in mod z^(K)
# For this particular loop iteration, we have sigma[l] and omega[l],
# and are computing sigma[K] and omega[K]
# First find Delta, the non-zero coefficient of z^(K) in
# (1 + s) * sigma[l]
# This delta is valid for l (this iteration) only
Delta = ( s2 * sigma[l] ).get_coefficient(l+1) # Delta is also known as the Discrepancy, and is always a scalar (not a polynomial).
# Make it a polynomial of degree 0, just for ease of computation with polynomials sigma and omega.
Delta = Polynomial(x0=Delta)
# Can now compute sigma[K] and omega[K] from
# sigma[l], omega[l], B[l], A[l], and Delta
sigma.append( sigma[l] - Delta * Z * B[l] )
omega.append( omega[l] - Delta * Z * A[l] )
# Now compute the next B and A
# There are two ways to do this
# This is based on a messy case analysis on the degrees of the four polynomials sigma, omega, A and B in order to minimize the degrees of A and B. For more infos, see https://www.cs.duke.edu/courses/spring10/cps296.3/decoding_rs_scribe.pdf
# In fact it ensures that the degree of the final polynomials aren't too large.
if Delta == ZERO or 2*L[l] > K \
or (2*L[l] == K and M[l] == 0):
# Rule A
B.append( Z * B[l] )
A.append( Z * A[l] )
L.append( L[l] )
M.append( M[l] )
elif (Delta != ZERO and 2*L[l] < K) \
or (2*L[l] == K and M[l] != 0):
# Rule B
B.append( sigma[l] // Delta )
A.append( omega[l] // Delta )
L.append( K - L[l] )
M.append( 1 - M[l] )
else:
raise Exception("Code shouldn't have gotten here")
return sigma[-1], omega[-1]
Polynomial 和 GF256int 分别是 2^8 上的多项式和 galois 域的通用实现。这些 类 经过单元测试，并且它们通常是错误证明。同样适用于 Reed-Solomon 的其余 encoding/decoding 方法，例如 Forney 和 Chien 搜索。可以在这里找到我在这里谈论的问题的完整代码和快速测试用例：http://codepad.org/l2Qi0y8o
这是一个示例输出：
Encoded message:
hello world�ꐙ�Ī`>
-------
Erasures decoding:
Erasure locator: 189x^5 + 88x^4 + 222x^3 + 33x^2 + 251x + 1
Syndrome: 149x^9 + 113x^8 + 29x^7 + 231x^6 + 210x^5 + 150x^4 + 192x^3 + 11x^2 + 41x
Sigma: 189x^5 + 88x^4 + 222x^3 + 33x^2 + 251x + 1
Symbols positions that were corrected: [19, 18, 17, 16, 15]
('Decoded message: ', 'hello world', '\xce\xea\x90\x99\x8d\xc4\xaa`>')
Correctly decoded: True
-------
Errors+Erasures decoding for the message with only erasures:
Erasure locator: 189x^5 + 88x^4 + 222x^3 + 33x^2 + 251x + 1
Syndrome: 149x^9 + 113x^8 + 29x^7 + 231x^6 + 210x^5 + 150x^4 + 192x^3 + 11x^2 + 41x
Sigma: 101x^10 + 139x^9 + 5x^8 + 14x^7 + 180x^6 + 148x^5 + 126x^4 + 135x^3 + 68x^2 + 155x + 1
Symbols positions that were corrected: [187, 141, 90, 19, 18, 17, 16, 15]
('Decoded message: ', '\xf4\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00.\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00P\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\xe3\xe6\xffO> world', '\xce\xea\x90\x99\x8d\xc4\xaa`>')
Correctly decoded: False
-------
Errors+Erasures decoding for the message with erasures and one error:
Erasure locator: 77x^4 + 96x^3 + 6x^2 + 206x + 1
Syndrome: 49x^9 + 107x^8 + x^7 + 109x^6 + 236x^5 + 15x^4 + 8x^3 + 133x^2 + 243x
Sigma: 38x^9 + 98x^8 + 239x^7 + 85x^6 + 32x^5 + 168x^4 + 92x^3 + 225x^2 + 22x + 1
Symbols positions that were corrected: [19, 18, 17, 16]
('Decoded message: ', "\xda\xe1'\xccA world", '\xce\xea\x90\x99\x8d\xc4\xaa`>')
Correctly decoded: False
这里，擦除解码总是正确的，因为它根本不使用 BM 来计算擦除定位器。通常，其他两个测试用例应该输出相同的 sigma，但它们根本就没有。
当你比较前两个测试用例时，问题来自 BM 的事实在这里是显而易见的：综合症和擦除定位器是相同的，但产生的西格玛完全不同（在第二个测试中，BM 是使用，而在第一个带有擦除的测试用例中，只有 BM 没有被调用）。
非常感谢您提供任何帮助或关于如何调试它的任何想法。请注意，您的答案可以是数学或代码，但请解释我的方法出了什么问题。
/EDIT: 仍然没有找到如何正确实现勘误表 BM 解码器的方法（请参阅下面我的回答）。赏金提供给任何可以解决问题（或至少指导我找到解决方案）的人。
/EDIT2: 愚蠢的我，抱歉，我刚刚重新阅读了架构，发现我错过了作业中的更改 L = r - L - erasures_count... 我有更新了代码以修复该问题并重新接受我的回答。
在阅读了大量的研究论文和书籍之后，我唯一找到答案的地方是在书中（readable online on Google Books，但没有 PDF 格式）：
"Algebraic codes for data transmission", Blahut, Richard E., 2003, Cambridge university press.
这里是本书的一些摘录，其中的细节是我实现的 Berlekamp-Massey 算法的准确描述（多项式运算的 matricial/vectorized 表示除外）：
这里是 Reed-Solomon 的勘误表（错误和擦除）Berlekamp-Massey 算法：
如您所见——与通常的描述相反，您只需要 初始化 Lambda，即错误定位器多项式，其值为先前计算的擦除定位器多项式 - - 您还需要跳过前 v 次迭代，其中 v 是擦除次数。请注意，这不等同于跳过最后 v 次迭代：您需要跳过前 v 次迭代，因为 r（迭代计数器，在我的实现中为 K）不仅用于计算迭代次数还要生成正确的差异系数 Delta。
这是经过修改以支持擦除以及最多 v+2*e <= (n-k):
的错误的结果代码
def _berlekamp_massey(self, s, k=None, erasures_loc=None, erasures_eval=None, erasures_count=0):
'''Computes and returns the errata (errors+erasures) locator polynomial (sigma) and the
error evaluator polynomial (omega) at the same time.
If the erasures locator is specified, we will return an errors-and-erasures locator polynomial and an errors-and-erasures evaluator polynomial, else it will compute only errors. With erasures in addition to errors, it can simultaneously decode up to v+2e <= (n-k) where v is the number of erasures and e the number of errors.
Mathematically speaking, this is equivalent to a spectral analysis (see Blahut, "Algebraic Codes for Data Transmission", 2003, chapter 7.6 Decoding in Time Domain).
The parameter s is the syndrome polynomial (syndromes encoded in a
generator function) as returned by _syndromes.
Notes:
The error polynomial:
E(x) = E_0 + E_1 x + ... + E_(n-1) x^(n-1)
j_1, j_2, ..., j_s are the error positions. (There are at most s
errors)
Error location X_i is defined: X_i = α^(j_i)
that is, the power of α (alpha) corresponding to the error location
Error magnitude Y_i is defined: E_(j_i)
that is, the coefficient in the error polynomial at position j_i
Error locator polynomial:
sigma(z) = Product( 1 - X_i * z, i=1..s )
roots are the reciprocals of the error locations
( 1/X_1, 1/X_2, ...)
Error evaluator polynomial omega(z) is here computed at the same time as sigma, but it can also be constructed afterwards using the syndrome and sigma (see _find_error_evaluator() method).
It can be seen that the algorithm tries to iteratively solve for the error locator polynomial by
solving one equation after another and updating the error locator polynomial. If it turns out that it
cannot solve the equation at some step, then it computes the error and weights it by the last
non-zero discriminant found, and delays the weighted result to increase the polynomial degree
by 1. Ref: "Reed Solomon Decoder: TMS320C64x Implementation" by Jagadeesh Sankaran, December 2000, Application Report SPRA686
The best paper I found describing the BM algorithm for errata (errors-and-erasures) evaluator computation is in "Algebraic Codes for Data Transmission", Richard E. Blahut, 2003.
'''
# For errors-and-erasures decoding, see: "Algebraic Codes for Data Transmission", Richard E. Blahut, 2003 and (but it's less complete): Blahut, Richard E. "Transform techniques for error control codes." IBM Journal of Research and development 23.3 (1979): 299-315. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.600&rep=rep1&type=pdf and also a MatLab implementation here: http://www.mathworks.com/matlabcentral/fileexchange/23567-reed-solomon-errors-and-erasures-decoder/content//RS_E_E_DEC.m
# also see: Blahut, Richard E. "A universal Reed-Solomon decoder." IBM Journal of Research and Development 28.2 (1984): 150-158. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.84.2084&rep=rep1&type=pdf
# and another good alternative book with concrete programming examples: Jiang, Yuan. A practical guide to error-control coding using Matlab. Artech House, 2010.
n = self.n
if not k: k = self.k
# Initialize, depending on if we include erasures or not:
if erasures_loc:
sigma = [ Polynomial(erasures_loc.coefficients) ] # copy erasures_loc by creating a new Polynomial, so that we initialize the errata locator polynomial with the erasures locator polynomial.
B = [ Polynomial(erasures_loc.coefficients) ]
omega = [ Polynomial(erasures_eval.coefficients) ] # to compute omega (the evaluator polynomial) at the same time, we also need to initialize it with the partial erasures evaluator polynomial
A = [ Polynomial(erasures_eval.coefficients) ] # TODO: fix the initial value of the evaluator support polynomial, because currently the final omega is not correct (it contains higher order terms that should be removed by the end of BM)
else:
sigma = [ Polynomial([GF256int(1)]) ] # error locator polynomial. Also called Lambda in other notations.
B = [ Polynomial([GF256int(1)]) ] # this is the error locator support/secondary polynomial, which is a funky way to say that it's just a temporary variable that will help us construct sigma, the error locator polynomial
omega = [ Polynomial([GF256int(1)]) ] # error evaluator polynomial. We don't need to initialize it with erasures_loc, it will still work, because Delta is computed using sigma, which itself is correctly initialized with erasures if needed.
A = [ Polynomial([GF256int(0)]) ] # this is the error evaluator support/secondary polynomial, to help us construct omega
L = [ 0 ] # update flag: necessary variable to check when updating is necessary and to check bounds (to avoid wrongly eliminating the higher order terms). For more infos, see https://www.cs.duke.edu/courses/spring11/cps296.3/decoding_rs.pdf
M = [ 0 ] # optional variable to check bounds (so that we do not mistakenly overwrite the higher order terms). This is not necessary, it's only an additional safe check. For more infos, see the presentation decoding_rs.pdf by Andrew Brown in the doc folder.
# Fix the syndrome shifting: when computing the syndrome, some implementations may prepend a 0 coefficient for the lowest degree term (the constant). This is a case of syndrome shifting, thus the syndrome will be bigger than the number of ecc symbols (I don't know what purpose serves this shifting). If that's the case, then we need to account for the syndrome shifting when we use the syndrome such as inside BM, by skipping those prepended coefficients.
# Another way to detect the shifting is to detect the 0 coefficients: by definition, a syndrome does not contain any 0 coefficient (except if there are no errors/erasures, in this case they are all 0). This however doesn't work with the modified Forney syndrome (that we do not use in this lib but it may be implemented in the future), which set to 0 the coefficients corresponding to erasures, leaving only the coefficients corresponding to errors.
synd_shift = 0
if len(s) > (n-k): synd_shift = len(s) - (n-k)
# Polynomial constants:
ONE = Polynomial(z0=GF256int(1))
ZERO = Polynomial(z0=GF256int(0))
Z = Polynomial(z1=GF256int(1)) # used to shift polynomials, simply multiply your poly * Z to shift
# Precaching
s2 = ONE + s
# Iteratively compute the polynomials n-k-erasures_count times. The last ones will be correct (since the algorithm refines the error/errata locator polynomial iteratively depending on the discrepancy, which is kind of a difference-from-correctness measure).
for l in xrange(0, n-k-erasures_count): # skip the first erasures_count iterations because we already computed the partial errata locator polynomial (by initializing with the erasures locator polynomial)
K = erasures_count+l+synd_shift # skip the FIRST erasures_count iterations (not the last iterations, that's very important!)
# Goal for each iteration: Compute sigma[l+1] and omega[l+1] such that
# (1 + s)*sigma[l] == omega[l] in mod z^(K)
# For this particular loop iteration, we have sigma[l] and omega[l],
# and are computing sigma[l+1] and omega[l+1]
# First find Delta, the non-zero coefficient of z^(K) in
# (1 + s) * sigma[l]
# Note that adding 1 to the syndrome s is not really necessary, you can do as well without.
# This delta is valid for l (this iteration) only
Delta = ( s2 * sigma[l] ).get_coefficient(K) # Delta is also known as the Discrepancy, and is always a scalar (not a polynomial).
# Make it a polynomial of degree 0, just for ease of computation with polynomials sigma and omega.
Delta = Polynomial(x0=Delta)
# Can now compute sigma[l+1] and omega[l+1] from
# sigma[l], omega[l], B[l], A[l], and Delta
sigma.append( sigma[l] - Delta * Z * B[l] )
omega.append( omega[l] - Delta * Z * A[l] )
# Now compute the next support polynomials B and A
# There are two ways to do this
# This is based on a messy case analysis on the degrees of the four polynomials sigma, omega, A and B in order to minimize the degrees of A and B. For more infos, see https://www.cs.duke.edu/courses/spring10/cps296.3/decoding_rs_scribe.pdf
# In fact it ensures that the degree of the final polynomials aren't too large.
if Delta == ZERO or 2*L[l] > K+erasures_count \
or (2*L[l] == K+erasures_count and M[l] == 0):
#if Delta == ZERO or len(sigma[l+1]) <= len(sigma[l]): # another way to compute when to update, and it doesn't require to maintain the update flag L
# Rule A
B.append( Z * B[l] )
A.append( Z * A[l] )
L.append( L[l] )
M.append( M[l] )
elif (Delta != ZERO and 2*L[l] < K+erasures_count) \
or (2*L[l] == K+erasures_count and M[l] != 0):
# elif Delta != ZERO and len(sigma[l+1]) > len(sigma[l]): # another way to compute when to update, and it doesn't require to maintain the update flag L
# Rule B
B.append( sigma[l] // Delta )
A.append( omega[l] // Delta )
L.append( K - L[l] ) # the update flag L is tricky: in Blahut's schema, it's mandatory to use `L = K - L - erasures_count` (and indeed in a previous draft of this function, if you forgot to do `- erasures_count` it would lead to correcting only 2*(errors+erasures) <= (n-k) instead of 2*errors+erasures <= (n-k)), but in this latest draft, this will lead to a wrong decoding in some cases where it should correctly decode! Thus you should try with and without `- erasures_count` to update L on your own implementation and see which one works OK without producing wrong decoding failures.
M.append( 1 - M[l] )
else:
raise Exception("Code shouldn't have gotten here")
# Hack to fix the simultaneous computation of omega, the errata evaluator polynomial: because A (the errata evaluator support polynomial) is not correctly initialized (I could not find any info in academic papers). So at the end, we get the correct errata evaluator polynomial omega + some higher order terms that should not be present, but since we know that sigma is always correct and the maximum degree should be the same as omega, we can fix omega by truncating too high order terms.
if omega[-1].degree > sigma[-1].degree: omega[-1] = Polynomial(omega[-1].coefficients[-(sigma[-1].degree+1):])
# Return the last result of the iterations (since BM compute iteratively, the last iteration being correct - it may already be before, but we're not sure)
return sigma[-1], omega[-1]
def _find_erasures_locator(self, erasures_pos):
'''Compute the erasures locator polynomial from the erasures positions (the positions must be relative to the x coefficient, eg: "hello worldxxxxxxxxx" is tampered to "h_ll_ worldxxxxxxxxx" with xxxxxxxxx being the ecc of length n-k=9, here the string positions are [1, 4], but the coefficients are reversed since the ecc characters are placed as the first coefficients of the polynomial, thus the coefficients of the erased characters are n-1 - [1, 4] = [18, 15] = erasures_loc to be specified as an argument.'''
# See: http://ocw.usu.edu/Electrical_and_Computer_Engineering/Error_Control_Coding/lecture7.pdf and Blahut, Richard E. "Transform techniques for error control codes." IBM Journal of Research and development 23.3 (1979): 299-315. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.600&rep=rep1&type=pdf and also a MatLab implementation here: http://www.mathworks.com/matlabcentral/fileexchange/23567-reed-solomon-errors-and-erasures-decoder/content//RS_E_E_DEC.m
erasures_loc = Polynomial([GF256int(1)]) # just to init because we will multiply, so it must be 1 so that the multiplication starts correctly without nulling any term
# erasures_loc is very simple to compute: erasures_loc = prod(1 - x*alpha[j]**i) for i in erasures_pos and where alpha is the alpha chosen to evaluate polynomials (here in this library it's gf(3)). To generate c*x where c is a constant, we simply generate a Polynomial([c, 0]) where 0 is the constant and c is positionned to be the coefficient for x^1. See https://en.wikipedia.org/wiki/Forney_algorithm#Erasures
for i in erasures_pos:
erasures_loc = erasures_loc * (Polynomial([GF256int(1)]) - Polynomial([GF256int(self.generator)**i, 0]))
return erasures_loc
注意：Sigma、Omega、A、B、L 和 M 都是多项式列表（因此我们保留了每次迭代计算的所有中间多项式的完整历史记录） .这当然可以优化，因为我们只需要 Sigma[l]、Sigma[l-1]、Omega[l]、Omega[l-1]、A[l]、B[l]、L[l] 和 M[l]（所以只有 Sigma 和 Omega 需要在内存中保留之前的迭代，其他变量不需要）。
注2：更新标志L是棘手的：在某些实现中，像在Blahut的模式中那样做会导致解码时错误的失败。在我过去的实现中，必须使用 L = K - L - erasures_count 来正确解码错误和擦除直到 Singleton 边界，但在我最新的实现中，我不得不使用 L = K - L（即使有擦除）以避免错误的解码失败。您应该在自己的实现中尝试这两种方法，看看哪一种不会产生任何错误的解码失败。有关更多信息，请参阅下面的问题。
此算法的唯一问题是它没有描述如何同时计算 Omega，即错误评估多项式（书中描述了如何仅针对错误初始化 Omega，而不是在解码错误和擦除时）。我尝试了几种变体和上述方法，但并不完全：最后，Omega 将包含本应取消的高阶项。可能是 Omega 或 A 错误评估器支持多项式，未使用合适的值进行初始化。
但是，您可以通过修剪太高阶项的 Omega 多项式来解决这个问题（因为它的次数应该与 Lambda/Sigma 相同）：
if omega[-1].degree > sigma[-1].degree: omega[-1] = Polynomial(omega[-1].coefficients[-(sigma[-1].degree+1):])
或者您完全可以在 BM 之后使用勘误表定位器 Lambda/Sigma 从头开始计算 Omega，它始终是正确计算的：
def _find_error_evaluator(self, synd, sigma, k=None):
'''Compute the error (or erasures if you supply sigma=erasures locator polynomial) evaluator polynomial Omega from the syndrome and the error/erasures/errata locator Sigma. Omega is already computed at the same time as Sigma inside the Berlekamp-Massey implemented above, but in case you modify Sigma, you can recompute Omega afterwards using this method, or just ensure that Omega computed by BM is correct given Sigma (as long as syndrome and sigma are correct, omega will be correct).'''
n = self.n
if not k: k = self.k
# Omega(x) = [ Synd(x) * Error_loc(x) ] mod x^(n-k+1) -- From Blahut, Algebraic codes for data transmission, 2003
return (synd * sigma) % Polynomial([GF256int(1)] + [GF256int(0)] * (n-k+1)) # Note that you should NOT do (1+Synd(x)) as can be seen in some books because this won't work with all primitive generators.
我正在 following question on CSTheory 中寻找更好的解决方案。
/编辑： 我将描述我遇到的一些问题以及如何解决它们：
- 不要忘记用擦除定位器多项式初始化错误定位器多项式（您可以轻松地从综合症和擦除位置计算）。
- 如果你只能解码错误并且只能完美擦除，但仅限于
2*errors + erasures <= (n-k)/2，那么你忘记跳过前 v 次迭代。
- 如果您可以解码擦除和错误但最多
2*(errors+erasures) <= (n-k)，那么您忘记更新 L 的分配：L = i+1 - L - erasures_count 而不是 L = i+1 - L。但在某些情况下，这实际上可能会使您的解码器失败，具体取决于您实现解码器的方式，请参阅下一点。
- 我的第一个解码器仅限于一个 generator/prime polynomial/fcr，但是当我将其更新为通用并添加严格的单元测试时，解码器在不应该的时候失败了。似乎 Blahut 上面的模式关于 L（更新标志）是错误的：它必须使用
L = K - L 而不是 L = K - L - erasures_count 进行更新，因为这有时会导致解码器失败，即使我们处于单例绑定下（因此我们应该正确解码！）。计算 L = K - L 不仅会解决这些解码问题，而且还会给出与不使用更新标志 L 的替代更新方式完全相同的结果（即条件 if Delta == ZERO or len(sigma[l+1]) <= len(sigma[l]):).但这很奇怪：在我过去的实现中，L = K - L - erasures_count 对于错误和擦除解码是强制性的，但现在它似乎会产生错误的失败。因此，您应该尝试使用和不使用您自己的实现，看看其中一个是否会为您带来错误的失败。
- 请注意，条件
2*L[l] > K 更改为 2*L[l] > K+erasures_count。如果一开始不添加条件 +erasures_count，您可能不会注意到任何副作用，但在某些情况下，解码会在不应该的情况下失败。
- 如果您只能修复一个错误或擦除，请检查您的条件是
2*L[l] > K+erasures_count 而不是 2*L[l] >= K+erasures_count（注意 > 而不是 >=）。
- 如果你可以纠正
2*errors + erasures <= (n-k-2)（刚好低于限制，例如，如果你有 10 个 ecc 符号，你只能纠正 4 个错误，而不是通常的 5 个错误）然后检查你的综合症和 BM 内的循环算法：如果综合症以常数项 x^0 的 0 系数开始（有时在书中建议），那么您的综合症会发生变化，然后 BM 内的循环必须从 1 开始并在 n-k+1 而不是 0:(n-k) 如果没有移位。
- 如果您可以更正除最后一个符号（最后一个 ecc 符号）之外的每个符号，那么请检查您的范围，尤其是在您的 Chien 搜索中：您不应该评估从 alpha^0 到 alpha^255 的错误定位器多项式但是从 alpha^1 到 alpha^256.
我参考了您的 python 代码并由 C++ 重写。
有效，您的信息和示例代码非常有用。
而且我发现错误的失败可能是由 M 值引起的。
According to "Algebraic codes for data transmission", the M value should not be a member of if-else case.
删除 M 后我没有遇到任何错误的失败。（或者只是还没有失败）
非常感谢您的知识分享。
// calculate C
Ref<ModulusPoly> T = C;
// M just for shift x
ArrayRef<int> m_temp(2);
m_temp[0]=1;
m_poly = new ModulusPoly(field_, m_temp);
// C = C - d*B*x
ArrayRef<int> d_temp(1);
d_temp[0] = d;
Ref<ModulusPoly> d_poly (new ModulusPoly(field_, d_temp));
d_poly = d_poly->multiply(m_poly);
d_poly = d_poly->multiply(B);
C = C->subtract(d_poly);
if(2*L<=n+e_size-1 && d!=0)
{
// b = d^-1
ArrayRef<int> b_temp(1);
b_temp[0] = field_.inverse(d);
b_poly = new ModulusPoly(field_, b_temp);
L = n-L-e_size;
// B = B*b = B*d^-1
B = T->multiply(b_poly);
}
else
{
// B = B*x
B = B->multiply(m_poly);
}
此答案是为了回应 gaborous 的评论而提供的。它没有显示如何修改 Berlekamp Massey 以处理擦除。相反，它显示了生成修改的 (Forney) 校正子的替代方案，该校正子消除了校正子中的擦除，之后修改的校正子可以与任何标准错误解码器算法一起使用以生成错误定位器多项式。然后将擦除和错误定位器多项式组合（相乘）以创建勘误定位器多项式。
这种方法不是最优的，因为有一些方法可以增强通用解码器以处理擦除和错误，但更通用。
用于在 C 中生成修改后的校正子的示例遗留代码。此代码中的 "vectors" 包括一个大小和一个数组。 vErsf 是一个与数据（代码字）大小相同的数组，非擦除位置为 0，擦除位置为 1。擦除被转换为擦除定位器多项式 (Root2Poly)，然后用于将标准校正子转换为修改后的 (Forney) 校正子。
typedef unsigned char ELEM; /* element type */
typedef struct{ /* vector structure */
ELEM size;
ELEM data[255];
}VECTOR;
static VECTOR vData; /* data */
static VECTOR vErsf; /* erasure flags (same size as data) */
static VECTOR vSyndromes; /* syndromes */
static VECTOR vErsLocators; /* erasure locators */
static VECTOR pErasures; /* erasure poly */
static VECTOR vModSyndromes; /* modified syndromes */
/*----------------------------------------------------------------------*/
/* Root2Poly(pPDst, pVSrc) convert roots into polynomial */
/*----------------------------------------------------------------------*/
static void Root2Poly(VECTOR *pPDst, VECTOR *pVSrc)
{
int i, j;
pPDst->size = pVSrc->size+1;
pPDst->data[0] = 1;
memset(&pPDst->data[1], 0, pVSrc->size*sizeof(ELEM));
for(j = 0; j < pVSrc->size; j++){
for(i = j; i >= 0; i--){
pPDst->data[i+1] = GFSub(pPDst->data[i+1],
GFMpy(pPDst->data[i], pVSrc->data[j]));}}
}
/*----------------------------------------------------------------------*/
/* GenErasures generate vErsLocat...and pErasures */
/* */
/* Scan vErsf right to left; for each non-zero flag byte, */
/* set vErsLocators to Alpha**power of corresponding location. */
/* Then convert these locators into a polynomial. */
/*----------------------------------------------------------------------*/
static int GenErasures(void)
{
int i, j;
ELEM eLcr; /* locator */
j = 0; /* j = index into vErs... */
eLcr = 1; /* init locator */
for(i = vErsf.size; i;){
i--;
if(vErsf.data[i]){ /* if byte flagged */
vErsLocators.data[j] = eLcr; /* set locator */
j++;
if(j > vGenRoots.size){ /* exit if too many */
return(1);}}
eLcr = GFMpy(eLcr, eAlpha);} /* bump locator */
vErsLocators.size = j; /* set size */
Root2Poly(&pErasures, &vErsLocators); /* generate poly */
return(0);
}
/*----------------------------------------------------------------------*/
/* GenModSyndromes Generate vModSyndromes */
/* */
/* generate modified syndromes */
/*----------------------------------------------------------------------*/
static void GenModSyndromes(void)
{
int iM; /* vModSyn.. index */
int iS; /* vSyn.. index */
int iP; /* pErs.. index */
vModSyndromes.size = vSyndromes.size-vErsLocators.size;
if(vErsLocators.size == 0){ /* if no erasures, copy */
memcpy(vModSyndromes.data, vSyndromes.data,
vSyndromes.size*sizeof(ELEM));
return;}
iS = 0;
for(iM = 0; iM < vModSyndromes.size; iM++){ /* modify */
vModSyndromes.data[iM] = 0;
iP = pErasures.size;
while(iP){
iP--;
vModSyndromes.data[iM] = GFAdd(vModSyndromes.data[iM],
GFMpy(vSyndromes.data[iS], pErasures.data[iP]));
iS++;}
iS -= pErasures.size-1;}
}
我正在尝试在 Python 中实现支持擦除和错误解码的 Reed-Solomon 编码器-解码器，这让我抓狂。
该实现目前支持仅解码错误或仅解码擦除，但不能同时解码（即使它低于 2*errors+erasures <= (n-k) 的理论界限）。
从 Blahut 的论文 (here and here) 看来，我们只需要用擦除定位器多项式初始化错误定位器多项式就可以隐式计算 Berlekamp-Massey 中的勘误定位器多项式。
这种方法对我部分有效：当我有 2*errors+erasures < (n-k)/2 时它有效，但实际上调试后它只有效，因为 BM 计算了一个错误定位器多项式，它得到了完全相同的值作为擦除定位器多项式（因为我们低于仅错误校正的限制），因此它通过伽罗瓦域被截断，我们最终得到擦除定位器多项式的正确值（至少我是这样理解的，我可能是错误的）。
然而，当我们超过 (n-k)/2 时，例如，如果 n = 20 和 k = 11，那么我们有 (n-k)=9 个我们可以纠正的擦除符号，如果我们输入 5 个擦除，那么 BM只是出错了。如果我们输入 4 个擦除 + 1 个错误（我们仍然远低于界限，因为我们有 2*错误+擦除 = 2+4 = 6 < 9），BM 仍然出错。
我实现的Berlekamp-Massey的具体算法可以在this presentation (pages 15-17), but a very similar description can be found here and here中找到，这里我附上一份数学描述：
现在，我几乎可以将这个数学算法精确地复制到 Python 代码中。我想要的是扩展它以支持擦除，我尝试通过使用擦除定位器初始化错误定位器 sigma:
def _berlekamp_massey(self, s, k=None, erasures_loc=None):
'''Computes and returns the error locator polynomial (sigma) and the
error evaluator polynomial (omega).
If the erasures locator is specified, we will return an errors-and-erasures locator polynomial and an errors-and-erasures evaluator polynomial.
The parameter s is the syndrome polynomial (syndromes encoded in a
generator function) as returned by _syndromes. Don't be confused with
the other s = (n-k)/2
Notes:
The error polynomial:
E(x) = E_0 + E_1 x + ... + E_(n-1) x^(n-1)
j_1, j_2, ..., j_s are the error positions. (There are at most s
errors)
Error location X_i is defined: X_i = a^(j_i)
that is, the power of a corresponding to the error location
Error magnitude Y_i is defined: E_(j_i)
that is, the coefficient in the error polynomial at position j_i
Error locator polynomial:
sigma(z) = Product( 1 - X_i * z, i=1..s )
roots are the reciprocals of the error locations
( 1/X_1, 1/X_2, ...)
Error evaluator polynomial omega(z) is here computed at the same time as sigma, but it can also be constructed afterwards using the syndrome and sigma (see _find_error_evaluator() method).
'''
# For errors-and-erasures decoding, see: Blahut, Richard E. "Transform techniques for error control codes." IBM Journal of Research and development 23.3 (1979): 299-315. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.600&rep=rep1&type=pdf and also a MatLab implementation here: http://www.mathworks.com/matlabcentral/fileexchange/23567-reed-solomon-errors-and-erasures-decoder/content//RS_E_E_DEC.m
# also see: Blahut, Richard E. "A universal Reed-Solomon decoder." IBM Journal of Research and Development 28.2 (1984): 150-158. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.84.2084&rep=rep1&type=pdf
# or alternatively see the reference book by Blahut: Blahut, Richard E. Theory and practice of error control codes. Addison-Wesley, 1983.
# and another good alternative book with concrete programming examples: Jiang, Yuan. A practical guide to error-control coding using Matlab. Artech House, 2010.
n = self.n
if not k: k = self.k
# Initialize:
if erasures_loc:
sigma = [ Polynomial(erasures_loc.coefficients) ] # copy erasures_loc by creating a new Polynomial
B = [ Polynomial(erasures_loc.coefficients) ]
else:
sigma = [ Polynomial([GF256int(1)]) ] # error locator polynomial. Also called Lambda in other notations.
B = [ Polynomial([GF256int(1)]) ] # this is the error locator support/secondary polynomial, which is a funky way to say that it's just a temporary variable that will help us construct sigma, the error locator polynomial
omega = [ Polynomial([GF256int(1)]) ] # error evaluator polynomial. We don't need to initialize it with erasures_loc, it will still work, because Delta is computed using sigma, which itself is correctly initialized with erasures if needed.
A = [ Polynomial([GF256int(0)]) ] # this is the error evaluator support/secondary polynomial, to help us construct omega
L = [ 0 ] # necessary variable to check bounds (to avoid wrongly eliminating the higher order terms). For more infos, see https://www.cs.duke.edu/courses/spring11/cps296.3/decoding_rs.pdf
M = [ 0 ] # optional variable to check bounds (so that we do not mistakenly overwrite the higher order terms). This is not necessary, it's only an additional safe check. For more infos, see the presentation decoding_rs.pdf by Andrew Brown in the doc folder.
# Polynomial constants:
ONE = Polynomial(z0=GF256int(1))
ZERO = Polynomial(z0=GF256int(0))
Z = Polynomial(z1=GF256int(1)) # used to shift polynomials, simply multiply your poly * Z to shift
s2 = ONE + s
# Iteratively compute the polynomials 2s times. The last ones will be
# correct
for l in xrange(0, n-k):
K = l+1
# Goal for each iteration: Compute sigma[K] and omega[K] such that
# (1 + s)*sigma[l] == omega[l] in mod z^(K)
# For this particular loop iteration, we have sigma[l] and omega[l],
# and are computing sigma[K] and omega[K]
# First find Delta, the non-zero coefficient of z^(K) in
# (1 + s) * sigma[l]
# This delta is valid for l (this iteration) only
Delta = ( s2 * sigma[l] ).get_coefficient(l+1) # Delta is also known as the Discrepancy, and is always a scalar (not a polynomial).
# Make it a polynomial of degree 0, just for ease of computation with polynomials sigma and omega.
Delta = Polynomial(x0=Delta)
# Can now compute sigma[K] and omega[K] from
# sigma[l], omega[l], B[l], A[l], and Delta
sigma.append( sigma[l] - Delta * Z * B[l] )
omega.append( omega[l] - Delta * Z * A[l] )
# Now compute the next B and A
# There are two ways to do this
# This is based on a messy case analysis on the degrees of the four polynomials sigma, omega, A and B in order to minimize the degrees of A and B. For more infos, see https://www.cs.duke.edu/courses/spring10/cps296.3/decoding_rs_scribe.pdf
# In fact it ensures that the degree of the final polynomials aren't too large.
if Delta == ZERO or 2*L[l] > K \
or (2*L[l] == K and M[l] == 0):
# Rule A
B.append( Z * B[l] )
A.append( Z * A[l] )
L.append( L[l] )
M.append( M[l] )
elif (Delta != ZERO and 2*L[l] < K) \
or (2*L[l] == K and M[l] != 0):
# Rule B
B.append( sigma[l] // Delta )
A.append( omega[l] // Delta )
L.append( K - L[l] )
M.append( 1 - M[l] )
else:
raise Exception("Code shouldn't have gotten here")
return sigma[-1], omega[-1]
Polynomial 和 GF256int 分别是 2^8 上的多项式和 galois 域的通用实现。这些 类 经过单元测试，并且它们通常是错误证明。同样适用于 Reed-Solomon 的其余 encoding/decoding 方法，例如 Forney 和 Chien 搜索。可以在这里找到我在这里谈论的问题的完整代码和快速测试用例：http://codepad.org/l2Qi0y8o
这是一个示例输出：
Encoded message:
hello world�ꐙ�Ī`>
-------
Erasures decoding:
Erasure locator: 189x^5 + 88x^4 + 222x^3 + 33x^2 + 251x + 1
Syndrome: 149x^9 + 113x^8 + 29x^7 + 231x^6 + 210x^5 + 150x^4 + 192x^3 + 11x^2 + 41x
Sigma: 189x^5 + 88x^4 + 222x^3 + 33x^2 + 251x + 1
Symbols positions that were corrected: [19, 18, 17, 16, 15]
('Decoded message: ', 'hello world', '\xce\xea\x90\x99\x8d\xc4\xaa`>')
Correctly decoded: True
-------
Errors+Erasures decoding for the message with only erasures:
Erasure locator: 189x^5 + 88x^4 + 222x^3 + 33x^2 + 251x + 1
Syndrome: 149x^9 + 113x^8 + 29x^7 + 231x^6 + 210x^5 + 150x^4 + 192x^3 + 11x^2 + 41x
Sigma: 101x^10 + 139x^9 + 5x^8 + 14x^7 + 180x^6 + 148x^5 + 126x^4 + 135x^3 + 68x^2 + 155x + 1
Symbols positions that were corrected: [187, 141, 90, 19, 18, 17, 16, 15]
('Decoded message: ', '\xf4\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00.\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00P\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\xe3\xe6\xffO> world', '\xce\xea\x90\x99\x8d\xc4\xaa`>')
Correctly decoded: False
-------
Errors+Erasures decoding for the message with erasures and one error:
Erasure locator: 77x^4 + 96x^3 + 6x^2 + 206x + 1
Syndrome: 49x^9 + 107x^8 + x^7 + 109x^6 + 236x^5 + 15x^4 + 8x^3 + 133x^2 + 243x
Sigma: 38x^9 + 98x^8 + 239x^7 + 85x^6 + 32x^5 + 168x^4 + 92x^3 + 225x^2 + 22x + 1
Symbols positions that were corrected: [19, 18, 17, 16]
('Decoded message: ', "\xda\xe1'\xccA world", '\xce\xea\x90\x99\x8d\xc4\xaa`>')
Correctly decoded: False
这里，擦除解码总是正确的，因为它根本不使用 BM 来计算擦除定位器。通常，其他两个测试用例应该输出相同的 sigma，但它们根本就没有。
当你比较前两个测试用例时，问题来自 BM 的事实在这里是显而易见的：综合症和擦除定位器是相同的，但产生的西格玛完全不同（在第二个测试中，BM 是使用，而在第一个带有擦除的测试用例中，只有 BM 没有被调用）。
非常感谢您提供任何帮助或关于如何调试它的任何想法。请注意，您的答案可以是数学或代码，但请解释我的方法出了什么问题。
/EDIT: 仍然没有找到如何正确实现勘误表 BM 解码器的方法（请参阅下面我的回答）。赏金提供给任何可以解决问题（或至少指导我找到解决方案）的人。
/EDIT2: 愚蠢的我，抱歉，我刚刚重新阅读了架构，发现我错过了作业中的更改 L = r - L - erasures_count... 我有更新了代码以修复该问题并重新接受我的回答。
在阅读了大量的研究论文和书籍之后，我唯一找到答案的地方是在书中（readable online on Google Books，但没有 PDF 格式）：
"Algebraic codes for data transmission", Blahut, Richard E., 2003, Cambridge university press.
这里是本书的一些摘录，其中的细节是我实现的 Berlekamp-Massey 算法的准确描述（多项式运算的 matricial/vectorized 表示除外）：
这里是 Reed-Solomon 的勘误表（错误和擦除）Berlekamp-Massey 算法：
如您所见——与通常的描述相反，您只需要 初始化 Lambda，即错误定位器多项式，其值为先前计算的擦除定位器多项式 - - 您还需要跳过前 v 次迭代，其中 v 是擦除次数。请注意，这不等同于跳过最后 v 次迭代：您需要跳过前 v 次迭代，因为 r（迭代计数器，在我的实现中为 K）不仅用于计算迭代次数还要生成正确的差异系数 Delta。
这是经过修改以支持擦除以及最多 v+2*e <= (n-k):
def _berlekamp_massey(self, s, k=None, erasures_loc=None, erasures_eval=None, erasures_count=0):
'''Computes and returns the errata (errors+erasures) locator polynomial (sigma) and the
error evaluator polynomial (omega) at the same time.
If the erasures locator is specified, we will return an errors-and-erasures locator polynomial and an errors-and-erasures evaluator polynomial, else it will compute only errors. With erasures in addition to errors, it can simultaneously decode up to v+2e <= (n-k) where v is the number of erasures and e the number of errors.
Mathematically speaking, this is equivalent to a spectral analysis (see Blahut, "Algebraic Codes for Data Transmission", 2003, chapter 7.6 Decoding in Time Domain).
The parameter s is the syndrome polynomial (syndromes encoded in a
generator function) as returned by _syndromes.
Notes:
The error polynomial:
E(x) = E_0 + E_1 x + ... + E_(n-1) x^(n-1)
j_1, j_2, ..., j_s are the error positions. (There are at most s
errors)
Error location X_i is defined: X_i = α^(j_i)
that is, the power of α (alpha) corresponding to the error location
Error magnitude Y_i is defined: E_(j_i)
that is, the coefficient in the error polynomial at position j_i
Error locator polynomial:
sigma(z) = Product( 1 - X_i * z, i=1..s )
roots are the reciprocals of the error locations
( 1/X_1, 1/X_2, ...)
Error evaluator polynomial omega(z) is here computed at the same time as sigma, but it can also be constructed afterwards using the syndrome and sigma (see _find_error_evaluator() method).
It can be seen that the algorithm tries to iteratively solve for the error locator polynomial by
solving one equation after another and updating the error locator polynomial. If it turns out that it
cannot solve the equation at some step, then it computes the error and weights it by the last
non-zero discriminant found, and delays the weighted result to increase the polynomial degree
by 1. Ref: "Reed Solomon Decoder: TMS320C64x Implementation" by Jagadeesh Sankaran, December 2000, Application Report SPRA686
The best paper I found describing the BM algorithm for errata (errors-and-erasures) evaluator computation is in "Algebraic Codes for Data Transmission", Richard E. Blahut, 2003.
'''
# For errors-and-erasures decoding, see: "Algebraic Codes for Data Transmission", Richard E. Blahut, 2003 and (but it's less complete): Blahut, Richard E. "Transform techniques for error control codes." IBM Journal of Research and development 23.3 (1979): 299-315. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.600&rep=rep1&type=pdf and also a MatLab implementation here: http://www.mathworks.com/matlabcentral/fileexchange/23567-reed-solomon-errors-and-erasures-decoder/content//RS_E_E_DEC.m
# also see: Blahut, Richard E. "A universal Reed-Solomon decoder." IBM Journal of Research and Development 28.2 (1984): 150-158. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.84.2084&rep=rep1&type=pdf
# and another good alternative book with concrete programming examples: Jiang, Yuan. A practical guide to error-control coding using Matlab. Artech House, 2010.
n = self.n
if not k: k = self.k
# Initialize, depending on if we include erasures or not:
if erasures_loc:
sigma = [ Polynomial(erasures_loc.coefficients) ] # copy erasures_loc by creating a new Polynomial, so that we initialize the errata locator polynomial with the erasures locator polynomial.
B = [ Polynomial(erasures_loc.coefficients) ]
omega = [ Polynomial(erasures_eval.coefficients) ] # to compute omega (the evaluator polynomial) at the same time, we also need to initialize it with the partial erasures evaluator polynomial
A = [ Polynomial(erasures_eval.coefficients) ] # TODO: fix the initial value of the evaluator support polynomial, because currently the final omega is not correct (it contains higher order terms that should be removed by the end of BM)
else:
sigma = [ Polynomial([GF256int(1)]) ] # error locator polynomial. Also called Lambda in other notations.
B = [ Polynomial([GF256int(1)]) ] # this is the error locator support/secondary polynomial, which is a funky way to say that it's just a temporary variable that will help us construct sigma, the error locator polynomial
omega = [ Polynomial([GF256int(1)]) ] # error evaluator polynomial. We don't need to initialize it with erasures_loc, it will still work, because Delta is computed using sigma, which itself is correctly initialized with erasures if needed.
A = [ Polynomial([GF256int(0)]) ] # this is the error evaluator support/secondary polynomial, to help us construct omega
L = [ 0 ] # update flag: necessary variable to check when updating is necessary and to check bounds (to avoid wrongly eliminating the higher order terms). For more infos, see https://www.cs.duke.edu/courses/spring11/cps296.3/decoding_rs.pdf
M = [ 0 ] # optional variable to check bounds (so that we do not mistakenly overwrite the higher order terms). This is not necessary, it's only an additional safe check. For more infos, see the presentation decoding_rs.pdf by Andrew Brown in the doc folder.
# Fix the syndrome shifting: when computing the syndrome, some implementations may prepend a 0 coefficient for the lowest degree term (the constant). This is a case of syndrome shifting, thus the syndrome will be bigger than the number of ecc symbols (I don't know what purpose serves this shifting). If that's the case, then we need to account for the syndrome shifting when we use the syndrome such as inside BM, by skipping those prepended coefficients.
# Another way to detect the shifting is to detect the 0 coefficients: by definition, a syndrome does not contain any 0 coefficient (except if there are no errors/erasures, in this case they are all 0). This however doesn't work with the modified Forney syndrome (that we do not use in this lib but it may be implemented in the future), which set to 0 the coefficients corresponding to erasures, leaving only the coefficients corresponding to errors.
synd_shift = 0
if len(s) > (n-k): synd_shift = len(s) - (n-k)
# Polynomial constants:
ONE = Polynomial(z0=GF256int(1))
ZERO = Polynomial(z0=GF256int(0))
Z = Polynomial(z1=GF256int(1)) # used to shift polynomials, simply multiply your poly * Z to shift
# Precaching
s2 = ONE + s
# Iteratively compute the polynomials n-k-erasures_count times. The last ones will be correct (since the algorithm refines the error/errata locator polynomial iteratively depending on the discrepancy, which is kind of a difference-from-correctness measure).
for l in xrange(0, n-k-erasures_count): # skip the first erasures_count iterations because we already computed the partial errata locator polynomial (by initializing with the erasures locator polynomial)
K = erasures_count+l+synd_shift # skip the FIRST erasures_count iterations (not the last iterations, that's very important!)
# Goal for each iteration: Compute sigma[l+1] and omega[l+1] such that
# (1 + s)*sigma[l] == omega[l] in mod z^(K)
# For this particular loop iteration, we have sigma[l] and omega[l],
# and are computing sigma[l+1] and omega[l+1]
# First find Delta, the non-zero coefficient of z^(K) in
# (1 + s) * sigma[l]
# Note that adding 1 to the syndrome s is not really necessary, you can do as well without.
# This delta is valid for l (this iteration) only
Delta = ( s2 * sigma[l] ).get_coefficient(K) # Delta is also known as the Discrepancy, and is always a scalar (not a polynomial).
# Make it a polynomial of degree 0, just for ease of computation with polynomials sigma and omega.
Delta = Polynomial(x0=Delta)
# Can now compute sigma[l+1] and omega[l+1] from
# sigma[l], omega[l], B[l], A[l], and Delta
sigma.append( sigma[l] - Delta * Z * B[l] )
omega.append( omega[l] - Delta * Z * A[l] )
# Now compute the next support polynomials B and A
# There are two ways to do this
# This is based on a messy case analysis on the degrees of the four polynomials sigma, omega, A and B in order to minimize the degrees of A and B. For more infos, see https://www.cs.duke.edu/courses/spring10/cps296.3/decoding_rs_scribe.pdf
# In fact it ensures that the degree of the final polynomials aren't too large.
if Delta == ZERO or 2*L[l] > K+erasures_count \
or (2*L[l] == K+erasures_count and M[l] == 0):
#if Delta == ZERO or len(sigma[l+1]) <= len(sigma[l]): # another way to compute when to update, and it doesn't require to maintain the update flag L
# Rule A
B.append( Z * B[l] )
A.append( Z * A[l] )
L.append( L[l] )
M.append( M[l] )
elif (Delta != ZERO and 2*L[l] < K+erasures_count) \
or (2*L[l] == K+erasures_count and M[l] != 0):
# elif Delta != ZERO and len(sigma[l+1]) > len(sigma[l]): # another way to compute when to update, and it doesn't require to maintain the update flag L
# Rule B
B.append( sigma[l] // Delta )
A.append( omega[l] // Delta )
L.append( K - L[l] ) # the update flag L is tricky: in Blahut's schema, it's mandatory to use `L = K - L - erasures_count` (and indeed in a previous draft of this function, if you forgot to do `- erasures_count` it would lead to correcting only 2*(errors+erasures) <= (n-k) instead of 2*errors+erasures <= (n-k)), but in this latest draft, this will lead to a wrong decoding in some cases where it should correctly decode! Thus you should try with and without `- erasures_count` to update L on your own implementation and see which one works OK without producing wrong decoding failures.
M.append( 1 - M[l] )
else:
raise Exception("Code shouldn't have gotten here")
# Hack to fix the simultaneous computation of omega, the errata evaluator polynomial: because A (the errata evaluator support polynomial) is not correctly initialized (I could not find any info in academic papers). So at the end, we get the correct errata evaluator polynomial omega + some higher order terms that should not be present, but since we know that sigma is always correct and the maximum degree should be the same as omega, we can fix omega by truncating too high order terms.
if omega[-1].degree > sigma[-1].degree: omega[-1] = Polynomial(omega[-1].coefficients[-(sigma[-1].degree+1):])
# Return the last result of the iterations (since BM compute iteratively, the last iteration being correct - it may already be before, but we're not sure)
return sigma[-1], omega[-1]
def _find_erasures_locator(self, erasures_pos):
'''Compute the erasures locator polynomial from the erasures positions (the positions must be relative to the x coefficient, eg: "hello worldxxxxxxxxx" is tampered to "h_ll_ worldxxxxxxxxx" with xxxxxxxxx being the ecc of length n-k=9, here the string positions are [1, 4], but the coefficients are reversed since the ecc characters are placed as the first coefficients of the polynomial, thus the coefficients of the erased characters are n-1 - [1, 4] = [18, 15] = erasures_loc to be specified as an argument.'''
# See: http://ocw.usu.edu/Electrical_and_Computer_Engineering/Error_Control_Coding/lecture7.pdf and Blahut, Richard E. "Transform techniques for error control codes." IBM Journal of Research and development 23.3 (1979): 299-315. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.600&rep=rep1&type=pdf and also a MatLab implementation here: http://www.mathworks.com/matlabcentral/fileexchange/23567-reed-solomon-errors-and-erasures-decoder/content//RS_E_E_DEC.m
erasures_loc = Polynomial([GF256int(1)]) # just to init because we will multiply, so it must be 1 so that the multiplication starts correctly without nulling any term
# erasures_loc is very simple to compute: erasures_loc = prod(1 - x*alpha[j]**i) for i in erasures_pos and where alpha is the alpha chosen to evaluate polynomials (here in this library it's gf(3)). To generate c*x where c is a constant, we simply generate a Polynomial([c, 0]) where 0 is the constant and c is positionned to be the coefficient for x^1. See https://en.wikipedia.org/wiki/Forney_algorithm#Erasures
for i in erasures_pos:
erasures_loc = erasures_loc * (Polynomial([GF256int(1)]) - Polynomial([GF256int(self.generator)**i, 0]))
return erasures_loc
注意：Sigma、Omega、A、B、L 和 M 都是多项式列表（因此我们保留了每次迭代计算的所有中间多项式的完整历史记录） .这当然可以优化，因为我们只需要 Sigma[l]、Sigma[l-1]、Omega[l]、Omega[l-1]、A[l]、B[l]、L[l] 和 M[l]（所以只有 Sigma 和 Omega 需要在内存中保留之前的迭代，其他变量不需要）。
注2：更新标志L是棘手的：在某些实现中，像在Blahut的模式中那样做会导致解码时错误的失败。在我过去的实现中，必须使用 L = K - L - erasures_count 来正确解码错误和擦除直到 Singleton 边界，但在我最新的实现中，我不得不使用 L = K - L（即使有擦除）以避免错误的解码失败。您应该在自己的实现中尝试这两种方法，看看哪一种不会产生任何错误的解码失败。有关更多信息，请参阅下面的问题。
此算法的唯一问题是它没有描述如何同时计算 Omega，即错误评估多项式（书中描述了如何仅针对错误初始化 Omega，而不是在解码错误和擦除时）。我尝试了几种变体和上述方法，但并不完全：最后，Omega 将包含本应取消的高阶项。可能是 Omega 或 A 错误评估器支持多项式，未使用合适的值进行初始化。
但是，您可以通过修剪太高阶项的 Omega 多项式来解决这个问题（因为它的次数应该与 Lambda/Sigma 相同）：
if omega[-1].degree > sigma[-1].degree: omega[-1] = Polynomial(omega[-1].coefficients[-(sigma[-1].degree+1):])
或者您完全可以在 BM 之后使用勘误表定位器 Lambda/Sigma 从头开始计算 Omega，它始终是正确计算的：
def _find_error_evaluator(self, synd, sigma, k=None):
'''Compute the error (or erasures if you supply sigma=erasures locator polynomial) evaluator polynomial Omega from the syndrome and the error/erasures/errata locator Sigma. Omega is already computed at the same time as Sigma inside the Berlekamp-Massey implemented above, but in case you modify Sigma, you can recompute Omega afterwards using this method, or just ensure that Omega computed by BM is correct given Sigma (as long as syndrome and sigma are correct, omega will be correct).'''
n = self.n
if not k: k = self.k
# Omega(x) = [ Synd(x) * Error_loc(x) ] mod x^(n-k+1) -- From Blahut, Algebraic codes for data transmission, 2003
return (synd * sigma) % Polynomial([GF256int(1)] + [GF256int(0)] * (n-k+1)) # Note that you should NOT do (1+Synd(x)) as can be seen in some books because this won't work with all primitive generators.
我正在 following question on CSTheory 中寻找更好的解决方案。
/编辑： 我将描述我遇到的一些问题以及如何解决它们：
- 不要忘记用擦除定位器多项式初始化错误定位器多项式（您可以轻松地从综合症和擦除位置计算）。
- 如果你只能解码错误并且只能完美擦除，但仅限于
2*errors + erasures <= (n-k)/2，那么你忘记跳过前 v 次迭代。 - 如果您可以解码擦除和错误但最多
2*(errors+erasures) <= (n-k)，那么您忘记更新 L 的分配：L = i+1 - L - erasures_count而不是L = i+1 - L。但在某些情况下，这实际上可能会使您的解码器失败，具体取决于您实现解码器的方式，请参阅下一点。 - 我的第一个解码器仅限于一个 generator/prime polynomial/fcr，但是当我将其更新为通用并添加严格的单元测试时，解码器在不应该的时候失败了。似乎 Blahut 上面的模式关于 L（更新标志）是错误的：它必须使用
L = K - L而不是L = K - L - erasures_count进行更新，因为这有时会导致解码器失败，即使我们处于单例绑定下（因此我们应该正确解码！）。计算L = K - L不仅会解决这些解码问题，而且还会给出与不使用更新标志 L 的替代更新方式完全相同的结果（即条件if Delta == ZERO or len(sigma[l+1]) <= len(sigma[l]):).但这很奇怪：在我过去的实现中，L = K - L - erasures_count对于错误和擦除解码是强制性的，但现在它似乎会产生错误的失败。因此，您应该尝试使用和不使用您自己的实现，看看其中一个是否会为您带来错误的失败。 - 请注意，条件
2*L[l] > K更改为2*L[l] > K+erasures_count。如果一开始不添加条件+erasures_count，您可能不会注意到任何副作用，但在某些情况下，解码会在不应该的情况下失败。 - 如果您只能修复一个错误或擦除，请检查您的条件是
2*L[l] > K+erasures_count而不是2*L[l] >= K+erasures_count（注意>而不是>=）。 - 如果你可以纠正
2*errors + erasures <= (n-k-2)（刚好低于限制，例如，如果你有 10 个 ecc 符号，你只能纠正 4 个错误，而不是通常的 5 个错误）然后检查你的综合症和 BM 内的循环算法：如果综合症以常数项 x^0 的 0 系数开始（有时在书中建议），那么您的综合症会发生变化，然后 BM 内的循环必须从1开始并在n-k+1而不是0:(n-k)如果没有移位。 - 如果您可以更正除最后一个符号（最后一个 ecc 符号）之外的每个符号，那么请检查您的范围，尤其是在您的 Chien 搜索中：您不应该评估从 alpha^0 到 alpha^255 的错误定位器多项式但是从 alpha^1 到 alpha^256.
我参考了您的 python 代码并由 C++ 重写。
有效，您的信息和示例代码非常有用。
而且我发现错误的失败可能是由 M 值引起的。
According to "Algebraic codes for data transmission", the
Mvalue should not be a member ofif-elsecase.
删除 M 后我没有遇到任何错误的失败。（或者只是还没有失败）
非常感谢您的知识分享。
// calculate C
Ref<ModulusPoly> T = C;
// M just for shift x
ArrayRef<int> m_temp(2);
m_temp[0]=1;
m_poly = new ModulusPoly(field_, m_temp);
// C = C - d*B*x
ArrayRef<int> d_temp(1);
d_temp[0] = d;
Ref<ModulusPoly> d_poly (new ModulusPoly(field_, d_temp));
d_poly = d_poly->multiply(m_poly);
d_poly = d_poly->multiply(B);
C = C->subtract(d_poly);
if(2*L<=n+e_size-1 && d!=0)
{
// b = d^-1
ArrayRef<int> b_temp(1);
b_temp[0] = field_.inverse(d);
b_poly = new ModulusPoly(field_, b_temp);
L = n-L-e_size;
// B = B*b = B*d^-1
B = T->multiply(b_poly);
}
else
{
// B = B*x
B = B->multiply(m_poly);
}
此答案是为了回应 gaborous 的评论而提供的。它没有显示如何修改 Berlekamp Massey 以处理擦除。相反，它显示了生成修改的 (Forney) 校正子的替代方案，该校正子消除了校正子中的擦除，之后修改的校正子可以与任何标准错误解码器算法一起使用以生成错误定位器多项式。然后将擦除和错误定位器多项式组合（相乘）以创建勘误定位器多项式。
这种方法不是最优的，因为有一些方法可以增强通用解码器以处理擦除和错误，但更通用。
用于在 C 中生成修改后的校正子的示例遗留代码。此代码中的 "vectors" 包括一个大小和一个数组。 vErsf 是一个与数据（代码字）大小相同的数组，非擦除位置为 0，擦除位置为 1。擦除被转换为擦除定位器多项式 (Root2Poly)，然后用于将标准校正子转换为修改后的 (Forney) 校正子。
typedef unsigned char ELEM; /* element type */
typedef struct{ /* vector structure */
ELEM size;
ELEM data[255];
}VECTOR;
static VECTOR vData; /* data */
static VECTOR vErsf; /* erasure flags (same size as data) */
static VECTOR vSyndromes; /* syndromes */
static VECTOR vErsLocators; /* erasure locators */
static VECTOR pErasures; /* erasure poly */
static VECTOR vModSyndromes; /* modified syndromes */
/*----------------------------------------------------------------------*/
/* Root2Poly(pPDst, pVSrc) convert roots into polynomial */
/*----------------------------------------------------------------------*/
static void Root2Poly(VECTOR *pPDst, VECTOR *pVSrc)
{
int i, j;
pPDst->size = pVSrc->size+1;
pPDst->data[0] = 1;
memset(&pPDst->data[1], 0, pVSrc->size*sizeof(ELEM));
for(j = 0; j < pVSrc->size; j++){
for(i = j; i >= 0; i--){
pPDst->data[i+1] = GFSub(pPDst->data[i+1],
GFMpy(pPDst->data[i], pVSrc->data[j]));}}
}
/*----------------------------------------------------------------------*/
/* GenErasures generate vErsLocat...and pErasures */
/* */
/* Scan vErsf right to left; for each non-zero flag byte, */
/* set vErsLocators to Alpha**power of corresponding location. */
/* Then convert these locators into a polynomial. */
/*----------------------------------------------------------------------*/
static int GenErasures(void)
{
int i, j;
ELEM eLcr; /* locator */
j = 0; /* j = index into vErs... */
eLcr = 1; /* init locator */
for(i = vErsf.size; i;){
i--;
if(vErsf.data[i]){ /* if byte flagged */
vErsLocators.data[j] = eLcr; /* set locator */
j++;
if(j > vGenRoots.size){ /* exit if too many */
return(1);}}
eLcr = GFMpy(eLcr, eAlpha);} /* bump locator */
vErsLocators.size = j; /* set size */
Root2Poly(&pErasures, &vErsLocators); /* generate poly */
return(0);
}
/*----------------------------------------------------------------------*/
/* GenModSyndromes Generate vModSyndromes */
/* */
/* generate modified syndromes */
/*----------------------------------------------------------------------*/
static void GenModSyndromes(void)
{
int iM; /* vModSyn.. index */
int iS; /* vSyn.. index */
int iP; /* pErs.. index */
vModSyndromes.size = vSyndromes.size-vErsLocators.size;
if(vErsLocators.size == 0){ /* if no erasures, copy */
memcpy(vModSyndromes.data, vSyndromes.data,
vSyndromes.size*sizeof(ELEM));
return;}
iS = 0;
for(iM = 0; iM < vModSyndromes.size; iM++){ /* modify */
vModSyndromes.data[iM] = 0;
iP = pErasures.size;
while(iP){
iP--;
vModSyndromes.data[iM] = GFAdd(vModSyndromes.data[iM],
GFMpy(vSyndromes.data[iS], pErasures.data[iP]));
iS++;}
iS -= pErasures.size-1;}
}