python - Errata (消去+エラー) Reed-Solomonデコーディング用のBerlekamp-Massey

Question

消去とエラーの両方のデコードをサポートする Python で Reed-Solomon エンコーダー/デコーダーを実装しようとしていますが、それが私を夢中にさせています。

実装は現在、エラーのみまたは消去のみのデコードをサポートしていますが、同時に両方をサポートすることはありません (たとえ 2*errors+erasures <= (nk) の理論的境界を下回っていたとしても)。

Blahut の論文 ( hereおよびhere ) から、Berlekamp-Massey 内のエラーロケータ多項式を暗黙的に計算するには、エラー ロケータ多項式を消去ロケータ多項式で初期化するだけでよいようです。

このアプローチは私にとって部分的に機能します: 2*errors+erasures < (nk)/2 の場合は機能しますが、実際にはデバッグ後にのみ機能します。これは、BM が消去ロケーター多項式とまったく同じ値を取得するエラー ロケーター多項式を計算するためです。 (エラーのみの修正の制限を下回っているため)、ガロア フィールドを介して切り捨てられ、消去ロケータ多項式の正しい値になります (少なくともそれが私が理解している方法であり、間違っている可能性があります)。

ただし、(nk)/2 を超えると、たとえば n = 20 で k = 11 の場合、(nk)=9 個の消去されたシンボルがあり、修正できます。5 個の消去をフィードすると、BM はうまくいきません。4 個の消去 + 1 個のエラー (2*エラー + 消去 = 2+4 = 6 < 9 であるため、まだ境界をはるかに下回っています) をフィードした場合でも、BM は正しく動作しません。

私が実装した Berlekamp-Massey の正確なアルゴリズムは、このプレゼンテーション(15 ～ 17 ページ) にありますが、非常によく似た説明がこことここにあります。ここに、数学的な説明のコピーを添付します。

これで、この数学的アルゴリズムを Python コードにほぼ正確に再現できました。私が望むのは、消去ロケーターでエラーロケーターシグマを初期化することによって試みた消去をサポートするように拡張することです:

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 上の多項式とガロア体の一般的な実装です。これらのクラスは単体テスト済みで、通常はバグが証明されています。Forney や Chien 検索など、Reed-Solomon の残りのエンコード/デコード方法についても同様です。ここで話している問題の簡単なテスト ケースを含む完全なコードは、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 をまったく使用しないため、イレージャー デコードは常に正しく行われます。通常、他の 2 つのテスト ケースは同じシグマを出力するはずですが、そうではありません。

問題が BM に起因するという事実は、最初の 2 つのテスト ケースを比較すると明らかです。シンドロームと消去ロケータは同じですが、結果のシグマはまったく異なります (2 番目のテストでは BM が使用され、消去のみの最初のテスト ケース BM は呼び出されません)。

これをデバッグする方法についての助けやアイデアをありがとう。あなたの答えは数学でもコードでもかまいませんが、私のアプローチで何がうまくいかなかったのか説明してください.

/編集:エラータ BM デコーダーを正しく実装する方法がまだ見つかりませんでした (以下の私の回答を参照してください)。この報奨金は、問題を解決できる (または少なくとも私を解決策に導いてくれる) 人に提供されます。

/EDIT2:愚かな私、申し訳ありませんが、スキーマを読み直したところ、割り当ての変更を見逃していることがわかりましたL = r - L - erasures_count...コードを更新して修正し、回答を再承認しました。

Answer 1

たくさんの研究論文や本を読んだ後、答えを見つけた唯一の場所は本の中にあります（Google Booksでオンラインで読むことができますが、PDFとしては利用できません）：

「データ伝送のための代数コード」、ブラハット、リチャード E.、2003 年、ケンブリッジ大学出版局。

以下は、この本の抜粋です。詳細は、私が実装した Berlekamp-Massey アルゴリズムの正確な説明です (多項式演算の行列/ベクトル化された表現を除く)。

そして、リードソロモンの正誤表 (エラーと消去) の Berlekamp-Massey アルゴリズムは次のとおりです。

ご覧のとおり、以前に計算された消去ロケーター多項式の値で、エラー ロケーター多項式である Lambda を初期化するだけでよいという通常の説明とは対照的に、最初の v 反復もスキップする必要があります。ここで、v は消去回数。最後の v iterations をスキップすることと同等ではないことに注意してください:最初の v iterations をスキップする必要があります。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

注: シグマ、オメガ、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 は注意が必要です。一部の実装では、ブラハットのスキーマと同じように実行すると、デコード時に誤ったエラーが発生します。L = K - L - erasures_count以前の実装では、Singleton バウンドまでエラーと消去の両方を正しくデコードするために使用することが必須でしたが、最新の実装でL = K - Lは (消去がある場合でも) 間違ったデコードの失敗を避けるために使用する必要がありました。独自の実装で両方を試して、どちらが間違ったデコードの失敗を引き起こさないかを確認する必要があります。詳細については、以下の問題を参照してください。

このアルゴリズムの唯一の問題は、エラー評価多項式である Omega を同時に計算する方法について説明していないことです (この本では、Omega をエラーに対してのみ初期化する方法について説明していますが、エラーと消去をデコードする場合については説明していません)。いくつかのバリエーションを試してみましたが、完全には機能しませんでした。最終的に、オメガには、キャンセルされるべき高次の用語が含まれます。おそらく Omega または A エラー エバリュエーター サポート多項式は、適切な値で初期化されていません。

ただし、次数が高すぎる項のオメガ多項式をトリミングすることで、これを修正できます (ラムダ/シグマと同じ次数を持つ必要があるため)。

if omega[-1].degree > sigma[-1].degree: omega[-1] = Polynomial(omega[-1].coefficients[-(sigma[-1].degree+1):])

または、常に正しく計算されるエラッタ ロケーター Lambda/Sigma を使用して、BM の後でゼロからオメガを完全に計算できます。

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.

CSTheory に関する次の質問で、より良い解決策を探しています。

/EDIT：私が抱えていたいくつかの問題とその修正方法について説明します：

• エラー位置多項式を消去位置多項式で初期化することを忘れないでください (シンドロームと消去位置から簡単に計算できます)。
• エラーのみをデコードし、消去のみを完璧に行うことができますが、に限定されている2*errors + erasures <= (n-k)/2場合は、最初の v 反復をスキップするのを忘れています。
• 消去とエラーの両方をデコードできますが、 までである場合は、 の代わりに2*(errors+erasures) <= (n-k)L: の割り当てを更新するのを忘れています。ただし、デコーダの実装方法によっては、これにより実際にデコーダが失敗する場合があります。次のポイントを参照してください。L = i+1 - L - erasures_countL = i+1 - L
• 私の最初のデコーダーは 1 つの生成器/素数多項式/fcr に限定されていましたが、それをユニバーサルに更新し、厳密な単体テストを追加すると、デコーダーは本来あるべきでないときに失敗しました。上記の Blahut のスキーマは L (更新フラグ) について間違っているようです: これは、 L = K - Land notを使用して更新する必要があります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が、場合によっては、デコードが失敗してはならない場合があります。
• エラーまたは消去を 1 つだけ修正できる場合は、条件がであること2*L[l] > K+erasures_countとそうでないことを確認してください ( の代わりに に2*L[l] >= K+erasures_count注意してください)。>>=
• 修正できる場合2*errors + erasures <= (n-k-2)(制限のすぐ下。たとえば、10 個の ecc シンボルがある場合、通常は 5 個ではなく 4 個のエラーしか修正できません)、シンドロームと BM アルゴリズム内のループをチェックします: シンドロームが係数 0 で始まる場合定数項 x^0 の場合 (本で推奨されることもあります)、シンドロームがシフトされ、BM 内のループは、シフトされていない場合ではなく、1で開始および終了する必要があります。n-k+10:(n-k)
• 最後のシンボル (最後の ecc シンボル) 以外のすべてのシンボルを修正できる場合は、特に Chien 検索で範囲を確認してください。エラー位置多項式を alpha^0 から alpha^255 までではなく、alpha^1 から alpha^1 まで評価する必要があります。アルファ^256。

Answer 2

私はあなたの python コードを参照し、によって書き直しましたC++。

あなたの情報とサンプルコードは本当に役に立ちます。

そして、値が原因で間違った失敗が発生する可能性があることがわかりましたM。

「データ伝送の代数コード」によると、M値はケースのメンバーであってはなりませんif-else。

が削除された後、間違った失敗はありませんでし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);
    }