Reed-Solomon correction in GF(2^16), seems work.

rs16.py

# Modified from
# http://en.wikiversity.org/wiki/Reed–Solomon_codes_for_coders

gf_exp = [0] * 65536 * 2 # Create list of 65536 elements. In Python 2.6+, consider using bytearray
gf_log = [0] * 65536
gf_exp[0] = 1
x = 1
for i in range(1, 65535):
    x <<= 1
    if x & 0x10000:
        x ^= 0x01002d
        # http://webhome.cs.uvic.ca/~dmaslov/gf2%5E16mult.html
    gf_exp[i] = x
    gf_log[x] = i
for i in range(65535, 65536 * 2):
     gf_exp[i] = gf_exp[i - 65535]
gf_log[gf_exp[65535]] = 65535 # Set last missing elements in gf_log[]


def gf_mul(x, y):
     if x==0 or y==0:
        return 0
     return gf_exp[gf_log[x] + gf_log[y]]


def gf_div(x, y):
     if y==0:
        raise ZeroDivisionError("division by zero")
     if x==0:
        return 0
     return gf_exp[gf_log[x] + 65535 - gf_log[y]]
     
     
def gf_poly_scale(p, x):
    # TODO: use list comprehension
    r = [0] * len(p)
    for i in range(0, len(p)):
        r[i] = gf_mul(p[i], x)
    return r
    
    
def gf_poly_add(p, q):
    r = [0] * max(len(p), len(q))
    for i in range(0, len(p)):
        r[i + len(r) - len(p)] = p[i]
    for i in range(0, len(q)):
        r[i + len(r) - len(q)] ^= q[i]
    return r
    
    
def gf_poly_mul(p, q):
    r = [0] * (len(p) + len(q) - 1)
    for j in range(0, len(q)):
        for i in range(0, len(p)):
            r[i + j] ^= gf_mul(p[i], q[j])
    return r
    
    
def gf_poly_eval(p,x):
    y = p[0]
    for i in range(1, len(p)):
        y = gf_mul(y,x) ^ p[i]
    return y
    
    
def rs_generator_poly(nsym):
    g = [1]
    for i in range(0,nsym):
        g = gf_poly_mul(g, [1, gf_exp[i]])
    return g
    
    
def rs_encode_msg(msg_in, nsym):
    gen = rs_generator_poly(nsym)
    msg_out = [0] * (len(msg_in) + nsym)
    for i in range(0, len(msg_in)):
        msg_out[i] = msg_in[i]
    for i in range(0, len(msg_in)):
        coef = msg_out[i]
        if coef != 0:
            for j in range(0, len(gen)):
                msg_out[i+j] ^= gf_mul(gen[j], coef)
    for i in range(0, len(msg_in)):
        msg_out[i] = msg_in[i]
    return msg_out
   
   
def rs_calc_syndromes(msg, nsym):
    synd = [0] * nsym
    for i in range(0, nsym):
        synd[i] = gf_poly_eval(msg, gf_exp[i])
    return synd
    
    
def rs_correct_errata(msg, synd, pos):
    # calculate error locator polynomial
    q = [1]
    for i in range(0, len(pos)):
        x = gf_exp[len(msg) - 1 - pos[i]]
        q = gf_poly_mul(q, [x, 1])
    # calculate error evaluator polynomial
    p = synd[0:len(pos)]
    p.reverse()
    p = gf_poly_mul(p, q)
    p = p[len(p)-len(pos):len(p)]
    # formal derivative of error locator eliminates even terms
    qprime = q[len(q)&1:len(q):2]
    # compute corrections
    for i in range(0, len(pos)):
        x = gf_exp[pos[i]+65536-len(msg)]
        y = gf_poly_eval(p, x)
        z = gf_poly_eval(qprime, gf_mul(x,x))
        msg[pos[i]] ^= gf_div(y, gf_mul(x,z))

test.py

#!/usr/bin/env python3
import random

from rs16 import *


def encode(msg, nsym):
    return rs_encode_msg(msg, nsym)

    
def decode(msg, pos, nsym):
    synd = rs_calc_syndromes(msg, nsym)
    out = msg[:]
    rs_correct_errata(out, synd, pos)
    return out[:-nsym]


msg = [500, 600, 700, 800, 900]
code = encode(msg, 3)
code[0], code[2] = [0, 0]
assert decode(code, [0, 2], 3) == msg

msg = [random.randrange(0, 0xffff) for i in range(1000)]
code = encode(msg, 300)
code[100:200] = [0] * 100
code[600:800] = [0] * 200
pos = list(range(100, 200)) + list(range(600, 800))
assert decode(code, pos, 300) == msg