implementation of Reed-Solomon Coding using a case in https://static.chunichi.co.jp/chunichi/pages/feature/science/galois_field_in_auto_factory.html

reed_solomon.sage

# https://static.chunichi.co.jp/chunichi/pages/feature/QR/galois_field_in_auto_factory.html
X = GF(2).polynomial_ring().gen()
poly = X ** 8 + X ** 4 + X ** 3 + X ** 2 + 1
F = GF(2 ** 8, name="a", modulus=poly)
R.<x> = PolynomialRing(F)


def tobin(x, n):
    x = Integer(x)
    nbits = x.nbits()
    assert nbits <= n
    return x.bits() + [0] * (n - nbits)

def frombin(v):
    return int("".join(map(str, v)), 2)

def toF(x):
    x = frombin(tobin(x, 8)[::-1])
    return F.fetch_int(x)

def fromF(x):
    x = x.integer_representation()
    x = frombin(tobin(x, 8)[::-1])
    return x


# f(x) = a(x)x**22 - r(x) = k(x) * g(x) = k(x) * (x-1)(x-2)(x-4)...(x-2**21)
data = [65, 36, 134, 247, 114, 5, 21, 34, 6, 54, 246, 70, 82, 7, 118, 247, 38, 183, 50, 224, 236, 17, 0]
PR.<x> = PolynomialRing(ZZ)
ax = 0
for d in data:
    ax *= x
    ax += toF(d)
gx = 1
for i in range(22):
    gx *= x - toF(2) ** i
rx = ax * x**22 % gx
fx_ = ax * x**22 - rx
for i in range(22):
    assert fx_(toF(2)**i) == toF(0)

read_data = [fromF(fx_[i]) for i in range(45)[::-1]]
# data[1], data[2] are incorrect
read_data[1] = 44
read_data[2] = 6
fx = 0
for d in read_data:
    fx *= x
    fx += toF(d)
Sx = 0
for i in range(22):
    Sx += fx(toF(2)**i) * x**i

f = x**22
g = Sx
qs = []
h = 0
while True:
    q, r = f.quo_rem(g)
    f, g = g, r
    qs.append(q)
    h += 1
    if r.degree() <= 11:
        break
f, g = 1, 0
for i in range(h):
    f, g = - qs[i] * f + g, f
gamma = toF(1) / f[0]
sigma = gamma * f
eta = (-1) ** h * gamma * r

for i in range(45):
    if sigma(toF(2)**-i) == toF(0):
        d_sigma = 0
        d_sigma = -toF(2)**-i * sigma / (x - toF(2)**-i)
        e = fromF(eta(toF(2)**-i) / -d_sigma(toF(2)**-i))
        print(i, e, f"{read_data[44-i]} -> {fromF(toF(read_data[44-i]) - toF(e))}")