docs: def _berlekamp_massey_fast(self, s, k=None, erasures_loc=None, erasures_eval=None, erasures_count=0):
          '''Faster implementation of errata (errors-and-erasures) Berlekamp-Massey.
          Returns the error locator polynomial (sigma) and the
          error evaluator polynomial (omega) with a faster implementation.
          '''
          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.
              sigmaprev = Polynomial(sigma.coefficients)
              B = Polynomial(sigma.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
              omegaprev = Polynomial(omega.coefficients)
              A = Polynomial(omega.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 = sigmaprev = Polynomial([GF2int(1)]) # error locator polynomial. Also called Lambda in other notations.
              sigmaprev = Polynomial([GF2int(1)]) # we need the previous iteration to compute the next value of the support polynomials
              B = Polynomial([GF2int(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 = omegaprev = Polynomial([GF2int(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.
              omegaprev = Polynomial([GF2int(1)])
              A = Polynomial([GF2int(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([GF2int(1)])
          ZERO = GF2int(0)
          Z = Polynomial([GF2int(1), GF2int(0)]) # 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 _range(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.mul_at(sigma, K) # Delta is also known as the Discrepancy, and is always a scalar (not a polynomial). We just need one coefficient at a specific degree, so we can optimize by computing only the polynomial multiplication at this term, and skip the others.
  
              # Can now compute sigma[l+1] and omega[l+1] from
              # sigma[l], omega[l], B[l], A[l], and Delta
              sigmaprev = sigma
              omegaprev = omega
              sigma = sigma - (Z * B).scale(Delta)
              omega = omega - (Z * A).scale(Delta)
  
              # 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 > K+erasures_count:
                  #or (2*L == K+erasures_count and M == 0):
              #if Delta == ZERO or len(sigma) <= len(sigmaprev): # another way to compute when to update, and it doesn't require to maintain the update flag L
                  # Rule A
                  B = Z * B
                  A = Z * A
                  #L = L
                  #M = M
  
              else:
              #elif (Delta != ZERO and 2*L < K+erasures_count) \
              #    or (2*L == K+erasures_count and M != 0):
              # elif Delta != ZERO and len(sigma) > len(sigmaprev): # another way to compute when to update, and it doesn't require to maintain the update flag L
                  # Rule B
                  B = sigmaprev.scale(Delta.inverse())
                  A = omegaprev.scale(Delta.inverse())
                  L = K - 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 = 1 - M
  
              #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.degree > sigma.degree: omega = Polynomial(omega.coefficients[-(sigma.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, omega