# encoding: UTF-8
from operator import mul
class RSCoder(object):
def __init__(self, n, k):
'''Creates a new Reed-Solomon Encoder/Decoder object configured with
the given n and k values.
n is the length of a codeword, must be less than 256
k is the length of the message, must be less than n
The code will have error correcting power s where 2s = n - k
The typical RSCoder is RSCoder(255, 223)
'''
if n < 0 or k < 0:
raise ValueError("n and k must be positive")
if not n < 256:
raise ValueError("n must be at most 255")
if not k < n:
raise ValueError("Codeword length n must be greater than message"
" length k")
self.n = n
self.k = k
self.g = {}
# Generate the generator polynomial for RS codes
# g(x) = (x-a^1)(x-a^2)...(x-a^(n-k))
# a is 3, a generator for GF(2^8)
g = Polynomial([GF256int(1)])
for alpha in xrange(0, 2): self.g[alpha] = g
for alpha in xrange(1,n+1):
p = Polynomial([GF256int(1), GF256int(3)**alpha])
g = g * p
self.g[n-alpha] = g
def encode(self, message, poly=False, k=None):
'''Encode a given string with reed-solomon encoding. Returns a byte
string with the k message bytes and n-k parity bytes at the end.
If a message is < k bytes long, it is assumed to be padded at the front
with null bytes.
The sequence returned is always n bytes long.
If poly is not False, returns the encoded Polynomial object instead of
the polynomial translated back to a string (useful for debugging)
'''
n = self.n
if not k: k = self.k
if len(message)>k:
raise ValueError("Message length is max %d. Message was %d" % (k,
len(message)))
# Encode message as a polynomial:
m = Polynomial([GF256int(ord(x)) for x in message])
# Shift polynomial up by n-k by multiplying by x^(n-k) to reserve the first n-k coefficients for the ecc. This effectively pad the message with \0 bytes for the lower coefficients (where the ecc will be placed).
mprime = m * Polynomial([GF256int(1)] + [GF256int(0)]*(n-k))
# mprime = q*g + b for some q
# so let's find b:
b = mprime % self.g[k]
# Subtract out b, so now c = q*g
c = mprime - b
# Since c is a multiple of g, it has (at least) n-k roots: a^1 through
# a^(n-k)
if not poly:
# Turn the polynomial c back into a byte string
return ''.join([chr(x) for x in c.coefficients]).rjust(n, "\0") # rjust is useful for the nostrip feature
else: return c
def check(self, r, k=None):
'''Check if there's any error in a message+ecc. Can be used before decoding, in addition to hashes to detect if the message was tampered, or after decoding to check that the message was fully recovered.
returns True/False
'''
n = self.n
if not k: k = self.k
g = self.g[k]
# Turn r into a polynomial
r = Polynomial([GF256int(ord(x)) for x in r])
# Compute the syndromes:
sz = self._syndromes(r, k=k)
# Checking that the syndrome is all 0 is sufficient to check if there are no more any errors in the decoded message
return sz.coefficients.count(GF256int(0)) == len(sz) # Faster than all()
def decode(self, r, nostrip=False, k=None, erasures_pos=None, only_erasures=False):
'''Given a received string or byte array r, attempts to decode it. If
it's a valid codeword, or if there are no more than (n-k)/2 errors, the
message is returned.
A message always has k bytes, if a message contained less it is left
padded with null bytes. When decoded, these leading null bytes are
stripped, but that can cause problems if decoding binary data. When
nostrip is True, messages returned are always k bytes long. This is
useful to make sure no data is lost when decoding binary data.
Theoretically, we have R(x) = C(x) + E(x) + V(x), where R is the received message, C is the correct message without errors nor erasures, E are the errors and V the erasures. Thus the goal is to compute E and V from R, so that we can compute: C(x) = R(x) - E(x) - V(x), and then we have our original message! The main problem of decoding is to solve the so-called Key Equation, here we use Berlekamp-Massey.
'''
n = self.n
if not k: k = self.k
# Turn r into a polynomial
rp = Polynomial([GF256int(ord(x)) for x in r])
# Compute the syndromes:
sz = self._syndromes(rp, k=k)
if sz.coefficients.count(GF256int(0)) == len(sz): # the code is already valid, there's nothing to do
# The last n-k bytes are parity
if nostrip:
return r[:-(n-k)], r[-(n-k):]
else:
return r[:-(n-k)].lstrip("\0"), r[-(n-k):]
# Erasures locator polynomial computation
erasures_loc = None
if erasures_pos:
# Convert string positions to coefficients positions for the algebra to work (see _find_erasures_locator(), ecc characters represent the first coefficients while the message is put last, so it's exactly the reverse of the string positions where the message is first and the ecc is last, thus it's just like if you read the message+ecc string in reverse)
erasures_pos = [len(r)-1-x for x in erasures_pos]
# Compute the erasure locator polynomial
erasures_loc = self._find_erasures_locator(erasures_pos)
if only_erasures:
sigma = erasures_loc
omega = self._find_error_evaluator(sz, sigma, k=k)
else:
# Find the error locator polynomial and error evaluator polynomial
# using the Berlekamp-Massey algorithm
sigma, omega = self._berlekamp_massey(sz, k=k, erasures_loc=erasures_loc)
#omega = self._find_error_evaluator(sz, sigma, k=k) # you can always double-check that computation of omega by BM is correct by analytically computing omega from sigma and syndrome.
# Now use Chien's procedure to find the error locations
# j is an array of integers representing the positions of the errors, 0
# being the rightmost byte
# X is a corresponding array of GF(2^8) values where X_i = alpha^(j_i)
X, j = self._chien_search(sigma)
# And finally, find the error magnitudes with Forney's Formula
# Y is an array of GF(2^8) values corresponding to the error magnitude
# at the position given by the j array
Y = self._forney(omega, X, k=k)
# Put the error and locations together to form the error polynomial
Elist = []
if len(Y) >= len(j):
for i in xrange(255):
if i in j:
Elist.append(Y[j.index(i)])
else:
Elist.append(GF256int(0))
E = Polynomial( Elist[::-1] )
else:
E = Polynomial()
# And we get our real codeword!
c = rp - E # Remember what we wrote above: R(x) = C(x) + E(x), so here to get back the original codeword C(x) = R(x) - E(x) ! (V(x) the erasures are here is included inside E(x))
#if len(c) > len(r): c = rp # failsafe: in case the correction went totally wrong (we repaired padded null bytes instead of the message! thus we end up with a longer message than what we should have), then we just return the uncorrected message. Note: we compare the length of c with r on purpose, that's not an error: if we compare with rp, if the first few characters were erased (null bytes) in r, then in rp the Polynomial will automatically skip them, thus the length will always be smaller in that case.
# Form it back into a string and return all but the last n-k bytes
ret = ''.join(chr(x) for x in c.coefficients[:-(n-k)])
ecc = ''.join(chr(x) for x in c.coefficients[-(n-k):]) # also return the corrected ecc, so that user can check()
# :-(
if nostrip:
# Polynomial objects don't store leading 0 coefficients, so we
# actually need to pad this to k bytes
return ret.rjust(k, "\0"), ecc
else:
return ret, ecc, sz, sigma, omega, j, X, erasures_loc
def _syndromes(self, r, k=None):
'''Given the received codeword r in the form of a Polynomial object,
computes the syndromes and returns the syndrome polynomial.
Mathematically, it's essentially equivalent to a Fourrier Transform (Chien search being the inverse).
'''
n = self.n
if not k: k = self.k
return Polynomial( [r.evaluate( GF256int(3)**l ) for l in xrange(n-k, 0, -1)] + [GF256int(0)] )
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**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.
for i in erasures_pos:
erasures_loc = erasures_loc * (Polynomial([GF256int(1)]) - Polynomial([GF256int(3)**i, 0]))
return erasures_loc
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]
def _find_error_evaluator(self, synd, sigma, k=None):
'''Compute the error evaluator polynomial Omega from the syndrome and the error 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. Can also be used to compute erasure evaluator polynomial if supplied with erasure locator polynomial instead of error locator poly.'''
n = self.n
if not k: k = self.k
# Omega(x) = [ (1 + Synd(x)) * Error_loc(x) ] mod x^(n-k)
ONE = Polynomial([GF256int(1)])
return ((ONE + synd) * sigma) % Polynomial([GF256int(1)] + [GF256int(0)] * (n-k))
def _chien_search(self, sigma):
'''Recall the definition of sigma, it has s roots. To find them, this
function evaluates sigma at all 255 non-zero points to find the roots
The inverse of the roots are X_i, the error locations
Returns a list X of error locations, and a corresponding list j of
error positions (the discrete log of the corresponding X value) The
lists are up to s elements large.
This is essentially an inverse Fourrier transform.
Important technical math note: This implementation is not actually
Chien's search. Chien's search is a way to evaluate the polynomial
such that each evaluation only takes constant time. This here simply
does 255 evaluations straight up, which is much less efficient.
Said differently, we simply do a bruteforce search by trial substitution to find the zeros of this polynomial, which identifies the error locations.
'''
# TODO: find a more efficient algorithm, this is the slowest part of the whole decoding process (~2.5 ms, while any other part is only ~400microsec). Could try the Pruned FFT from "Simple Algorithms for BCH Decoding", by Jonathan Hong and Martin Vetterli, IEEE Transactions on Communications, Vol.43, No.8, August 1995
X = []
j = []
p = GF256int(3)
# Try for each possible location
for l in xrange(1,256):
# These evaluations could be more efficient, but oh well
if sigma.evaluate( p**l ) == 0: # If it's 0, then bingo! It's an error location
X.append( p**(-l) )
# This is different than the notes, I think the notes were in error
# Notes said j values were just l, when it's actually 255-l
j.append(255 - l)
return X, j
def _forney(self, omega, X, k=None):
'''Computes the error magnitudes. Works also with erasures and errors+erasures beyond the (n-k)//2 bound, here the bound is 2*e+v <= (n-k-1) with e the number of errors and v the number of erasures.'''
# XXX Is floor division okay here? Should this be ceiling?
if not k: k = self.k
Y = [] # the final result, the error/erasures polynomial (contain the values that we should minus on the received message to get the repaired message)
for l, Xl in enumerate(X):
# Compute the sequence product
Xl_inv = Xl.inverse()
Xlength = len(X)
prod = [1 - Xl_inv * X[j] for j in xrange(Xlength) if j != l]
prod = reduce(mul, prod, 1)
# Compute Yl
# This is a more faithful translation of the theoretical equation contrary to the old forney method. Here it is exactly copy/pasted from the included presentation decoding_rs.pdf: Yl = omega(Xl.inverse()) / prod(1 - Xj*Xl.inverse()) for j in len(X) (in the paper it's for j in s, but it's useless when len(X) < s because we compute neutral terms 1 for nothing, and wrong when correcting more than s erasures or erasures+errors since it prevents computing all required terms).
# Thus here this method works with erasures too because firstly we fixed the equation to be like the theoretical one (don't know why it was modified in _old_forney(), if it's an optimization, it doesn't enhance anything), and secondly because we removed the product bound on s, which prevented computing errors and erasures above the s=(n-k)//2 bound.
Yl = omega.evaluate(Xl_inv) / prod
Y.append(Yl)
return Y
#-------------------------------------------------------------------------
# POLYNOMIAL and FF (Galois Field) Classes, from here the implementation is generic and verified by various unit testing
#-------------------------------------------------------------------------
from cStringIO import StringIO
from itertools import izip
class Polynomial(object):
'''Completely general polynomial class.
Polynomial objects are immutable.
Implementation note: while this class is mostly agnostic to the type of
coefficients used (as long as they support the usual mathematical
operations), the Polynomial class still assumes the additive identity and
multiplicative identity are 0 and 1 respectively. If you're doing math over
some strange field or using non-numbers as coefficients, this class will
need to be modified.'''
def __init__(self, coefficients=None, **sparse):
'''
There are three ways to initialize a Polynomial object.
1) With a list, tuple, or other iterable, creates a polynomial using
the items as coefficients in order of decreasing power
2) With keyword arguments such as for example x3=5, sets the
coefficient of x^3 to be 5
3) With no arguments, creates an empty polynomial, equivalent to
Polynomial((0,))
>>> print Polynomial((5, 0, 0, 0, 0, 0))
5x^5
>>> print Polynomial(x32=5, x64=8)
8x^64 + 5x^32
>>> print Polynomial(x5=5, x9=4, x0=2)
4x^9 + 5x^5 + 2
'''
if coefficients is not None and sparse:
raise TypeError("Specify coefficients list /or/ keyword terms, not"
" both")
if coefficients is not None:
# Polynomial( [1, 2, 3, ...] )
#if isinstance(coefficients, tuple): coefficients = list(coefficients)
# Expunge any leading 0 coefficients
while len(coefficients) > 0 and coefficients[0] == 0:
coefficients.pop(0)
if not coefficients:
coefficients.append(0)
self.coefficients = coefficients
elif sparse:
# Polynomial(x32=...)
powers = sparse.keys()
powers.sort(reverse=1)
# Not catching possible exceptions from the following line, let
# them bubble up.
highest = int(powers[0][1:])
coefficients = [0] * (highest+1)
for power, coeff in sparse.iteritems():
power = int(power[1:])
coefficients[highest - power] = coeff
self.coefficients = coefficients
else:
# Polynomial()
self.coefficients = [0]
# In any case, compute the degree of the polynomial (=the maximum degree)
self.degree = len(self.coefficients)-1
def __len__(self):
"""Returns the number of terms in the polynomial"""
return self.degree+1
# return len(self.coefficients)
def get_degree(self, poly=None):
"""Returns the degree of the polynomial"""
if not poly:
return self.degree
#return len(self.coefficients) - 1
elif poly and hasattr("coefficients", poly):
return len(poly.coefficients) - 1
else:
while poly and poly[-1] == 0:
poly.pop() # normalize
return len(poly)-1
def __add__(self, other):
diff = len(self) - len(other)
t1 = [0] * (-diff) + self.coefficients
t2 = [0] * diff + other.coefficients
return self.__class__([x+y for x,y in izip(t1, t2)])
def __neg__(self):
if self[0].__class__.__name__ == "GF256int": # optimization: -GF256int(x) == GF256int(x), so it's useless to do a loop in this case
return self
else:
return self.__class__([-x for x in self.coefficients])
def __sub__(self, other):
return self + -other
def __mul__(self, other):
terms = [0] * (len(self) + len(other))
l1l2 = self.degree + other.degree
for i1, c1 in enumerate(self.coefficients):
if c1 == 0:
# Optimization
continue
for i2, c2 in enumerate(other.coefficients):
if c2 == 0:
continue
else:
terms[ -(l1l2-(i1+i2)+1) ] += c1*c2
return self.__class__(terms)
def mul_at(self, other, k):
'''Compute the multiplication between two polynomials only at the specified coefficient (this is a lot cheaper than doing the full polynomial multiplication and then extract only the required coefficient)'''
if k > (self.degree + other.degree): return 0 # optimization: if the required coefficient is above the maximum coefficient of the resulting polynomial, we can already predict that and just return 0
term = 0
for i in xrange(min(len(self), len(other))):
coef1 = self.coefficients[-(k-i+1)]
coef2 = other.coefficients[-(i+1)]
if coef1 == 0 or coef2 == 0: continue # optimization, skip if any coef is 0, the multiplication will bring nothing but 0
term += coef1 * coef2
return term
def scale(self, other):
'''Multiply a polynomial with a scalar'''
return self.__class__([self.coefficients[i] * other for i in xrange(len(self))])
def __floordiv__(self, other):
return divmod(self, other)[0]
def __mod__(self, other):
return divmod(self, other)[1]
def __divmod__(dividend, divisor):
'''Implementation of the Polynomial Long Division, without recursion. Polynomial Long Division is very similar to a simple division of integers, see purplemath.com. Implementation inspired by the pseudo-code from Rosettacode.org'''
'''Pseudocode:
degree(P):
return the index of the last non-zero element of P;
if all elements are 0, return -inf
polynomial_long_division(N, D) returns (q, r):
// N, D, q, r are vectors
if degree(D) < 0 then error
if degree(N) >= degree(D) then
q = 0
while degree(N) >= degree(D)
d = D shifted right by (degree(N) - degree(D))
q[degree(N) - degree(D)] = N(degree(N)) / d(degree(d))
// by construction, degree(d) = degree(N) of course
d = d * q[degree(N) - degree(D)]
N = N - d
endwhile
r = N
else
q = 0
r = N
endif
return (q, r)
'''
class_ = dividend.__class__
# See how many times the highest order term
# of the divisor can go into the highest order term of the dividend
dividend_power = dividend.degree
dividend_coefficient = dividend[0]
divisor_power = divisor.degree
divisor_coefficient = divisor[0]
if divisor_power < 0:
raise ZeroDivisionError
elif dividend_power < divisor_power:
# Doesn't divide at all, return 0 for the quotient and the entire
# dividend as the remainder
quotient = class_()
remainder = dividend
else: # dividend_power >= divisor_power
quotient = class_() # init the quotient array
# init the remainder to the dividend, and we will divide it sucessively by the quotient major coefficient
remainder = dividend
remainder_power = dividend_power
remainder_coefficient = dividend_coefficient
quotient_power = remainder_power - divisor_power # need to set at least 1 just to start the loop. Warning if set to remainder_power - divisor_power: because it may skip the loop altogether (and we want to at least do one iteration to set the quotient)
# Compute how many times the highest order term in the divisor goes into the dividend
while quotient_power >= 0 and remainder.coefficients != [0]: # Until there's no remainder left (or the remainder cannot be divided anymore by the divisor)
quotient_coefficient = remainder_coefficient / divisor_coefficient
q = class_( [quotient_coefficient] + [0] * quotient_power ) # construct an array with only the quotient major coefficient (we divide the remainder only with the major coeff)
quotient = quotient + q # add the coeff to the full quotient
remainder = remainder - q * divisor # divide the remainder with the major coeff quotient multiplied by the divisor, this gives us the new remainder
remainder_power = remainder.degree # compute the new remainder degree
remainder_coefficient = remainder[0] # Compute the new remainder coefficient
quotient_power = remainder_power - divisor_power
return quotient, remainder
def __eq__(self, other):
return self.coefficients == other.coefficients
def __ne__(self, other):
return self.coefficients != other.coefficients
def __hash__(self):
return hash(self.coefficients)
def __repr__(self):
n = self.__class__.__name__
return "%s(%r)" % (n, self.coefficients)
def __str__(self):
buf = StringIO()
l = len(self) - 1
for i, c in enumerate(self.coefficients):
if not c and i > 0:
continue
power = l - i
if c == 1 and power != 0:
c = ""
if power > 1:
buf.write("%sx^%s" % (c, power))
elif power == 1:
buf.write("%sx" % c)
else:
buf.write("%s" % c)
buf.write(" + ")
return buf.getvalue()[:-3]
def evaluate(self, x):
'''Evaluate this polynomial at value x, returning the result (which is the sum of all evaluations at each term).'''
# Faster implementation using Horner's Scheme
y = self[0]
for i in xrange(1, len(self)):
y = y * x + self.coefficients[i]
return y
def evaluate_array(self, x):
'''Simple way of evaluating a polynomial at value x, but here we return both the full array (evaluated at each polynomial position) and the sum'''
x_gf = self.coefficients[0].__class__(x)
arr = [self.coefficients[-i]*x_gf**(i-1) for i in xrange(len(self), 0, -1)]
return arr, sum(arr)
def derive(self):
'''Compute the formal derivative of the polynomial: sum(i*coeff[i] x^(i-1))'''
L = len(self)-1
return Polynomial( [(L-i) * self[i] for i in xrange(0, len(self)-1)] )
def get_coefficient(self, degree):
'''Returns the coefficient of the specified term'''
if degree > self.degree:
return 0
else:
return self.coefficients[-(degree+1)]
def __iter__(self):
return iter(self.coefficients)
def __getitem__(self, slice):
return self.coefficients[slice]
def __setitem__(self, key, item):
self.coefficients[key] = item
# Exponent table for 3, a generator for GF(256)
GF256int_exptable = [1, 3, 5, 15, 17, 51, 85, 255, 26, 46, 114, 150, 161, 248, 19,
53, 95, 225, 56, 72, 216, 115, 149, 164, 247, 2, 6, 10, 30, 34,
102, 170, 229, 52, 92, 228, 55, 89, 235, 38, 106, 190, 217, 112,
144, 171, 230, 49, 83, 245, 4, 12, 20, 60, 68, 204, 79, 209, 104,
184, 211, 110, 178, 205, 76, 212, 103, 169, 224, 59, 77, 215, 98,
166, 241, 8, 24, 40, 120, 136, 131, 158, 185, 208, 107, 189, 220,
127, 129, 152, 179, 206, 73, 219, 118, 154, 181, 196, 87, 249, 16,
48, 80, 240, 11, 29, 39, 105, 187, 214, 97, 163, 254, 25, 43, 125,
135, 146, 173, 236, 47, 113, 147, 174, 233, 32, 96, 160, 251, 22,
58, 78, 210, 109, 183, 194, 93, 231, 50, 86, 250, 21, 63, 65, 195,
94, 226, 61, 71, 201, 64, 192, 91, 237, 44, 116, 156, 191, 218,
117, 159, 186, 213, 100, 172, 239, 42, 126, 130, 157, 188, 223,
122, 142, 137, 128, 155, 182, 193, 88, 232, 35, 101, 175, 234, 37,
111, 177, 200, 67, 197, 84, 252, 31, 33, 99, 165, 244, 7, 9, 27,
45, 119, 153, 176, 203, 70, 202, 69, 207, 74, 222, 121, 139, 134,
145, 168, 227, 62, 66, 198, 81, 243, 14, 18, 54, 90, 238, 41, 123,
141, 140, 143, 138, 133, 148, 167, 242, 13, 23, 57, 75, 221, 124,
132, 151, 162, 253, 28, 36, 108, 180, 199, 82, 246, 1]
# Logarithm table, base 3
GF256int_logtable = [None, 0, 25, 1, 50, 2, 26, 198, 75, 199, 27, 104, 51, 238, 223,
3, 100, 4, 224, 14, 52, 141, 129, 239, 76, 113, 8, 200, 248, 105,
28, 193, 125, 194, 29, 181, 249, 185, 39, 106, 77, 228, 166, 114,
154, 201, 9, 120, 101, 47, 138, 5, 33, 15, 225, 36, 18, 240, 130,
69, 53, 147, 218, 142, 150, 143, 219, 189, 54, 208, 206, 148, 19,
92, 210, 241, 64, 70, 131, 56, 102, 221, 253, 48, 191, 6, 139, 98,
179, 37, 226, 152, 34, 136, 145, 16, 126, 110, 72, 195, 163, 182,
30, 66, 58, 107, 40, 84, 250, 133, 61, 186, 43, 121, 10, 21, 155,
159, 94, 202, 78, 212, 172, 229, 243, 115, 167, 87, 175, 88, 168,
80, 244, 234, 214, 116, 79, 174, 233, 213, 231, 230, 173, 232, 44,
215, 117, 122, 235, 22, 11, 245, 89, 203, 95, 176, 156, 169, 81,
160, 127, 12, 246, 111, 23, 196, 73, 236, 216, 67, 31, 45, 164,
118, 123, 183, 204, 187, 62, 90, 251, 96, 177, 134, 59, 82, 161,
108, 170, 85, 41, 157, 151, 178, 135, 144, 97, 190, 220, 252, 188,
149, 207, 205, 55, 63, 91, 209, 83, 57, 132, 60, 65, 162, 109, 71,
20, 42, 158, 93, 86, 242, 211, 171, 68, 17, 146, 217, 35, 32, 46,
137, 180, 124, 184, 38, 119, 153, 227, 165, 103, 74, 237, 222, 197,
49, 254, 24, 13, 99, 140, 128, 192, 247, 112, 7]
class GF256int(int):
'''Instances of this object are elements of the field GF(2^8)
Instances are integers in the range 0 to 255
This field is defined using the irreducable polynomial
x^8 + x^4 + x^3 + x + 1
and using 3 as the generator for the exponent table and log table.
'''
def __add__(a, b):
"Addition in GF(2^8) is the xor of the two"
return GF256int(a ^ b)
__sub__ = __add__
__radd__ = __add__
__rsub__ = __add__
def __neg__(self):
return self
def __mul__(a, b):
"Multiplication in GF(2^8)"
if a == 0 or b == 0:
return GF256int(0)
x = GF256int_logtable[a]
y = GF256int_logtable[b]
z = (x + y) % 255
return GF256int(GF256int_exptable[z])
__rmul__ = __mul__
def __pow__(self, power, modulo=None):
if isinstance(power, GF256int):
raise TypeError("Raising a Field element to another Field element is not defined. power must be a regular integer")
x = GF256int_logtable[self]
z = (x * power) % 255
return GF256int(GF256int_exptable[z])
def inverse(self):
e = GF256int_logtable[self]
return GF256int(GF256int_exptable[255 - e])
def __div__(self, other):
return self * GF256int(other).inverse()
def __rdiv__(self, other):
return self.inverse() * other
def __repr__(self):
n = self.__class__.__name__
return "%s(%r)" % (n, int(self))
# Simple test case
if __name__ == "__main__":
# Encoding hello world
rsman = RSCoder(20, 11)
mes = 'hello world'
mesecc = rsman.encode(mes)
print("Encoded message:")
print(mesecc)
# Erase some characters, beyond the errors limit
ii = 5
mesecc2 = "\x00" * ii + mesecc[ii:]
mesecc3 = "\x00" * (ii-1) + "A" + mesecc[ii:] # error one character, the rest is erased
erasures_pos = [i for i in xrange(len(mesecc2)) if mesecc2[i] == "\x00"]
# Correct only using erasures decoding
eras_r1, eras_r2, eras_synd, eras_sigma, eras_omega, eras_j, eras_X, eras_eras_loc = rsman.decode(mesecc2, erasures_pos=erasures_pos, only_erasures=True)
# Correct using errors+erasures decoding
r1, r2, synd, sigma, omega, j, X, eras_loc = rsman.decode(mesecc2, erasures_pos=erasures_pos)
# Correct using errors+erasures decoding the message containing erasures and one error
err_r1, err_r2, err_synd, err_sigma, err_omega, err_j, err_X, err_eras_loc = rsman.decode(mesecc3, erasures_pos=erasures_pos[:-1])
# Print results
print("-------")
print("Erasures decoding:")
print "Erasure locator:", eras_eras_loc
print "Syndrome:", eras_synd
print "Sigma:", eras_sigma
print "Symbols positions that were corrected:", eras_j
print("Decoded message: ", eras_r1, eras_r2)
print "Correctly decoded: ", rsman.check(eras_r1 + eras_r2)
print("-------")
print("Errors+Erasures decoding for the message with only erasures:")
print "Erasure locator:", eras_loc
print "Syndrome:", synd
print "Sigma:", sigma
print "Symbols positions that were corrected:", j
print("Decoded message: ", r1, r2)
print "Correctly decoded: ", rsman.check(r1 + r2)
print("-------")
print("Errors+Erasures decoding for the message with erasures and one error:")
print "Erasure locator:", err_eras_loc
print "Syndrome:", err_synd
print "Sigma:", err_sigma
print "Symbols positions that were corrected:", err_j
print("Decoded message: ", err_r1, err_r2)
print "Correctly decoded: ", rsman.check(err_r1 + err_r2)