def decode(self, r, nostrip=False, k=None, erasures_pos=None, only_erasures=False, return_string=True):
        '''Given a received string or byte array or list r of values between
        0 and gf2_charac, attempts to decode it. If it's a valid codeword, or
        if there are no more than (n-k)/2 errors, the repaired 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.
        
        When stated in the language of spectral estimation, consists of a Fourier transform (syndrome computer), followed by a spectral analysis (Berlekamp-Massey or Euclidian algorithm), followed by an inverse Fourier transform (Chien search).
        (see Blahut, "Algebraic Codes for Data Transmission", 2003, chapter 7.6 Decoding in Time Domain).
        '''
        n = self.n
        if not k: k = self.k

        # If we were given a string, convert to a list (important to support fields above 2^8)
        if isinstance(r, _str):
            r = [ord(x) for x in r]

        # Turn r into a polynomial
        rp = Polynomial([GF2int(x) for x in r])

        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]
            # Set erasures characters to null bytes
            # Note that you can just leave the original characters as they are, you don't need to set erased characters to null bytes for the decoding to work, but note that it won't help either (ie, fake erasures, meaning characters that were detected as erasures but actually aren't, will still "consume" one ecc symbol, even if you don't set them to null byte, this is because the syndrome is limited to n-k and thus you can't decode above this bound without a clever trick).
            # Example string containing a fake erasure: "hello sam" -> "ooooo sam" with erasures_pos = [0, 1, 2, 3, 4]. Here in fact the last erasure is fake because the original character also was "o" so if we detect "o" as an erasure, we will end up with one fake erasure. But setting it to null byte or not, it will still use up one ecc symbol, it will always be counted as a real erasure. If you're below the n-k bound, then the doceding will be ok. If you're above, then you can't do anything, the decoding won't work. Maybe todo: try to find a clever list decoding algorithm to account for fake erasures....
            # Note: commented out so that the resulting omega (error evaluator polynomial) is the same as the erasure evaluator polynomial when decoding the same number of errors or erasures (ie, decoding 3 erasures only will give the same result as 3 errors only, with of course the errors/erasures on the same characters).
            #for erasure in erasures_pos:
                #rp[erasure] = GF2int(0)

        # Compute the syndromes:
        sz = self._syndromes(rp, k=k)

        if sz.coefficients.count(GF2int(0)) == len(sz): # the code is already valid, there's nothing to do
            # The last n-k bytes are parity
            ret = r[:-(n-k)]
            ecc = r[-(n-k):]
            if not nostrip:
                ret = self._list_lstrip(r[:-(n-k)], 0)
            if return_string and self.gf2_charac < 256:
                ret = self._list2str(ret)
                ecc = self._list2str(ecc) 
            return ret, ecc

        # Erasures locator polynomial computation
        erasures_loc = None
        erasures_eval = None
        erasures_count = 0
        if erasures_pos:
            erasures_count = len(erasures_pos)
            # Compute the erasure locator polynomial
            erasures_loc = self._find_erasures_locator(erasures_pos)
            # Compute the erasure evaluator polynomial
            erasures_eval = self._find_error_evaluator(sz, erasures_loc, k=k)

        if only_erasures:
            sigma = erasures_loc
            omega = erasures_eval
        else:
            # Find the error locator polynomial and error evaluator polynomial
            # using the Berlekamp-Massey algorithm
            # if erasures were supplied, BM will generate the errata (errors-and-erasures) locator and evaluator polynomials
            sigma, omega = self._berlekamp_massey(sz, k=k, erasures_loc=erasures_loc, erasures_eval=erasures_eval, erasures_count=erasures_count)
            omega = self._find_error_evaluator(sz, sigma, k=k) # we want to make sure that omega is correct (we know that sigma is always correct, but omega not really)

        # 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)

        # Sanity check: Cannot guarantee correct decoding of more than n-k errata (Singleton Bound, n-k being the minimum distance), and we cannot even check if it's correct (the syndrome will always be all 0 if we try to decode above the bound), thus it's better to just return the input as-is.
        if len(j) > n-k:
            ret = r[:-(n-k)]
            ecc = r[-(n-k):]
            if not nostrip:
                ret = self._list_lstrip(r[:-(n-k)], 0)
            if return_string and self.gf2_charac < 256:
                ret = self._list2str(ret)
                ecc = self._list2str(ecc) 
            return ret, ecc

        # 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)

        # Put the error and locations together to form the error polynomial
        # Note that an alternative would be to compute the error-spectrum polynomial E(x) which satisfies E(x)*Sigma(x) = 0 (mod x^n - 1) = Omega(x)(x^n - 1) -- see Blahut, Algebraic codes for data transmission
        Elist = [GF2int(0)] * self.gf2_charac
        if len(Y) >= len(j): # failsafe: if the number of erratas is higher than the number of coefficients in the magnitude polynomial, we failed!
            for i in _range(self.gf2_charac): # FIXME? is this really necessary to go to self.gf2_charac? len(rp) wouldn't be just enough? (since the goal is anyway to substract E to rp)
                if i in j:
                    Elist[i] = Y[j.index(i)]
            E = Polynomial( Elist[::-1] ) # reverse the list because we used the coefficient degrees (j) instead of the error positions
        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.

        # Split the polynomial into two parts: the corrected message and the corrected ecc
        ret = c.coefficients[:-(n-k)]
        ecc = c.coefficients[-(n-k):]

        if nostrip:
            # Polynomial objects don't store leading 0 coefficients, so we
            # actually need to pad this to k bytes
            ret = self._list_rjust(ret, k, 0)

        if return_string and self.gf2_charac < 256: # automatically disable return_string if the field is above 255 (chr would fail, so it's up to the user to define the mapping)
            # Form it back into a string 
            ret = self._list2str(ret)
            ecc = self._list2str(ecc)

        return ret, ecc