ssdv/rs8.c

319 wiersze
9.1 KiB
C

/* General purpose Reed-Solomon codec for 8-bit symbols or less
* Copyright 2003 Phil Karn, KA9Q
* May be used under the terms of the GNU Lesser General Public License (LGPL)
*
* This version tweaked by Philip Heron <phil@sanslogic.co.uk>
*/
#include <string.h>
#include "rs8.h"
static const uint8_t ALPHA_TO[] = {
0x01,0x02,0x04,0x08,0x10,0x20,0x40,0x80,0x87,0x89,0x95,0xAD,0xDD,0x3D,0x7A,0xF4,
0x6F,0xDE,0x3B,0x76,0xEC,0x5F,0xBE,0xFB,0x71,0xE2,0x43,0x86,0x8B,0x91,0xA5,0xCD,
0x1D,0x3A,0x74,0xE8,0x57,0xAE,0xDB,0x31,0x62,0xC4,0x0F,0x1E,0x3C,0x78,0xF0,0x67,
0xCE,0x1B,0x36,0x6C,0xD8,0x37,0x6E,0xDC,0x3F,0x7E,0xFC,0x7F,0xFE,0x7B,0xF6,0x6B,
0xD6,0x2B,0x56,0xAC,0xDF,0x39,0x72,0xE4,0x4F,0x9E,0xBB,0xF1,0x65,0xCA,0x13,0x26,
0x4C,0x98,0xB7,0xE9,0x55,0xAA,0xD3,0x21,0x42,0x84,0x8F,0x99,0xB5,0xED,0x5D,0xBA,
0xF3,0x61,0xC2,0x03,0x06,0x0C,0x18,0x30,0x60,0xC0,0x07,0x0E,0x1C,0x38,0x70,0xE0,
0x47,0x8E,0x9B,0xB1,0xE5,0x4D,0x9A,0xB3,0xE1,0x45,0x8A,0x93,0xA1,0xC5,0x0D,0x1A,
0x34,0x68,0xD0,0x27,0x4E,0x9C,0xBF,0xF9,0x75,0xEA,0x53,0xA6,0xCB,0x11,0x22,0x44,
0x88,0x97,0xA9,0xD5,0x2D,0x5A,0xB4,0xEF,0x59,0xB2,0xE3,0x41,0x82,0x83,0x81,0x85,
0x8D,0x9D,0xBD,0xFD,0x7D,0xFA,0x73,0xE6,0x4B,0x96,0xAB,0xD1,0x25,0x4A,0x94,0xAF,
0xD9,0x35,0x6A,0xD4,0x2F,0x5E,0xBC,0xFF,0x79,0xF2,0x63,0xC6,0x0B,0x16,0x2C,0x58,
0xB0,0xE7,0x49,0x92,0xA3,0xC1,0x05,0x0A,0x14,0x28,0x50,0xA0,0xC7,0x09,0x12,0x24,
0x48,0x90,0xA7,0xC9,0x15,0x2A,0x54,0xA8,0xD7,0x29,0x52,0xA4,0xCF,0x19,0x32,0x64,
0xC8,0x17,0x2E,0x5C,0xB8,0xF7,0x69,0xD2,0x23,0x46,0x8C,0x9F,0xB9,0xF5,0x6D,0xDA,
0x33,0x66,0xCC,0x1F,0x3E,0x7C,0xF8,0x77,0xEE,0x5B,0xB6,0xEB,0x51,0xA2,0xC3,0x00,
};
static const uint8_t INDEX_OF[] = {
0xFF,0x00,0x01,0x63,0x02,0xC6,0x64,0x6A,0x03,0xCD,0xC7,0xBC,0x65,0x7E,0x6B,0x2A,
0x04,0x8D,0xCE,0x4E,0xC8,0xD4,0xBD,0xE1,0x66,0xDD,0x7F,0x31,0x6C,0x20,0x2B,0xF3,
0x05,0x57,0x8E,0xE8,0xCF,0xAC,0x4F,0x83,0xC9,0xD9,0xD5,0x41,0xBE,0x94,0xE2,0xB4,
0x67,0x27,0xDE,0xF0,0x80,0xB1,0x32,0x35,0x6D,0x45,0x21,0x12,0x2C,0x0D,0xF4,0x38,
0x06,0x9B,0x58,0x1A,0x8F,0x79,0xE9,0x70,0xD0,0xC2,0xAD,0xA8,0x50,0x75,0x84,0x48,
0xCA,0xFC,0xDA,0x8A,0xD6,0x54,0x42,0x24,0xBF,0x98,0x95,0xF9,0xE3,0x5E,0xB5,0x15,
0x68,0x61,0x28,0xBA,0xDF,0x4C,0xF1,0x2F,0x81,0xE6,0xB2,0x3F,0x33,0xEE,0x36,0x10,
0x6E,0x18,0x46,0xA6,0x22,0x88,0x13,0xF7,0x2D,0xB8,0x0E,0x3D,0xF5,0xA4,0x39,0x3B,
0x07,0x9E,0x9C,0x9D,0x59,0x9F,0x1B,0x08,0x90,0x09,0x7A,0x1C,0xEA,0xA0,0x71,0x5A,
0xD1,0x1D,0xC3,0x7B,0xAE,0x0A,0xA9,0x91,0x51,0x5B,0x76,0x72,0x85,0xA1,0x49,0xEB,
0xCB,0x7C,0xFD,0xC4,0xDB,0x1E,0x8B,0xD2,0xD7,0x92,0x55,0xAA,0x43,0x0B,0x25,0xAF,
0xC0,0x73,0x99,0x77,0x96,0x5C,0xFA,0x52,0xE4,0xEC,0x5F,0x4A,0xB6,0xA2,0x16,0x86,
0x69,0xC5,0x62,0xFE,0x29,0x7D,0xBB,0xCC,0xE0,0xD3,0x4D,0x8C,0xF2,0x1F,0x30,0xDC,
0x82,0xAB,0xE7,0x56,0xB3,0x93,0x40,0xD8,0x34,0xB0,0xEF,0x26,0x37,0x0C,0x11,0x44,
0x6F,0x78,0x19,0x9A,0x47,0x74,0xA7,0xC1,0x23,0x53,0x89,0xFB,0x14,0x5D,0xF8,0x97,
0x2E,0x4B,0xB9,0x60,0x0F,0xED,0x3E,0xE5,0xF6,0x87,0xA5,0x17,0x3A,0xA3,0x3C,0xB7,
};
static const uint8_t GENPOLY[] = {
0x00,0xF9,0x3B,0x42,0x04,0x2B,0x7E,0xFB,0x61,0x1E,0x03,0xD5,0x32,0x42,0xAA,0x05,
0x18,0x05,0xAA,0x42,0x32,0xD5,0x03,0x1E,0x61,0xFB,0x7E,0x2B,0x04,0x42,0x3B,0xF9,
0x00,
};
static inline int mod255(int x)
{
while(x >= 255)
{
x -= 255;
x = (x >> 8) + (x & 255);
}
return(x);
}
#define MODNN(x) mod255(x)
#define MM (8)
#define NN (255)
#define NROOTS (32)
#define FCR (112)
#define PRIM (11)
#define IPRIM (116)
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define A0 (NN) /* Special reserved value encoding zero in index form */
/* Portable C version */
void encode_rs_8(uint8_t *data, uint8_t *parity, int pad)
{
int i, j;
uint8_t feedback;
memset(parity, 0, NROOTS * sizeof(uint8_t));
for(i = 0; i < NN - NROOTS - pad; i++)
{
feedback = INDEX_OF[data[i] ^ parity[0]];
if(feedback != A0) /* feedback term is non-zero */
{
for(j = 1; j < NROOTS; j++)
parity[j] ^= ALPHA_TO[mod255(feedback + GENPOLY[NROOTS - j])];
}
/* Shift */
memmove(&parity[0], &parity[1], sizeof(uint8_t) * (NROOTS - 1));
if(feedback != A0)
parity[NROOTS - 1] = ALPHA_TO[mod255(feedback + GENPOLY[0])];
else
parity[NROOTS - 1] = 0;
}
}
int decode_rs_8(uint8_t *data, int *eras_pos, int no_eras, int pad)
{
int deg_lambda, el, deg_omega;
int i, j, r, k;
uint8_t u, q, tmp, num1, num2, den, discr_r;
uint8_t lambda[NROOTS + 1], s[NROOTS]; /* Err+Eras Locator poly
* and syndrome poly */
uint8_t b[NROOTS + 1], t[NROOTS + 1], omega[NROOTS + 1];
uint8_t root[NROOTS], reg[NROOTS + 1], loc[NROOTS];
int syn_error, count;
if(pad < 0 || pad > 222) return(-1);
/* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
for(i = 0; i < NROOTS; i++) s[i] = data[0];
for(j = 1; j < NN - pad; j++)
{
for(i = 0; i < NROOTS; i++)
{
if(s[i] == 0) s[i] = data[j];
else s[i] = data[j] ^ ALPHA_TO[MODNN(INDEX_OF[s[i]] + (FCR + i) * PRIM)];
}
}
/* Convert syndromes to index form, checking for nonzero condition */
syn_error = 0;
for(i = 0; i < NROOTS; i++)
{
syn_error |= s[i];
s[i] = INDEX_OF[s[i]];
}
if(!syn_error)
{
/* if syndrome is zero, data[] is a codeword and there are no
* errors to correct. So return data[] unmodified
*/
count = 0;
goto finish;
}
memset(&lambda[1], 0, NROOTS * sizeof(lambda[0]));
lambda[0] = 1;
if(no_eras > 0)
{
/* Init lambda to be the erasure locator polynomial */
lambda[1] = ALPHA_TO[MODNN(PRIM * (NN - 1 - eras_pos[0]))];
for(i = 1; i < no_eras; i++)
{
u = MODNN(PRIM * (NN - 1 - eras_pos[i]));
for(j = i + 1; j > 0; j--)
{
tmp = INDEX_OF[lambda[j - 1]];
if(tmp != A0) lambda[j] ^= ALPHA_TO[MODNN(u + tmp)];
}
}
}
for(i = 0; i < NROOTS + 1; i++)
b[i] = INDEX_OF[lambda[i]];
/*
* Begin Berlekamp-Massey algorithm to determine error+erasure
* locator polynomial
*/
r = no_eras;
el = no_eras;
while(++r <= NROOTS) /* r is the step number */
{
/* Compute discrepancy at the r-th step in poly-form */
discr_r = 0;
for(i = 0; i < r; i++)
{
if((lambda[i] != 0) && (s[r - i - 1] != A0))
{
discr_r ^= ALPHA_TO[MODNN(INDEX_OF[lambda[i]] + s[r - i - 1])];
}
}
discr_r = INDEX_OF[discr_r]; /* Index form */
if(discr_r == A0)
{
/* 2 lines below: B(x) <-- x*B(x) */
memmove(&b[1], b, NROOTS * sizeof(b[0]));
b[0] = A0;
}
else
{
/* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
t[0] = lambda[0];
for(i = 0; i < NROOTS; i++)
{
if(b[i] != A0)
t[i + 1] = lambda[i + 1] ^ ALPHA_TO[MODNN(discr_r + b[i])];
else
t[i + 1] = lambda[i + 1];
}
if(2 * el <= r + no_eras - 1)
{
el = r + no_eras - el;
/*
* 2 lines below: B(x) <-- inv(discr_r) *
* lambda(x)
*/
for(i = 0; i <= NROOTS; i++)
b[i] = (lambda[i] == 0) ? A0 : MODNN(INDEX_OF[lambda[i]] - discr_r + NN);
}
else
{
/* 2 lines below: B(x) <-- x*B(x) */
memmove(&b[1], b, NROOTS * sizeof(b[0]));
b[0] = A0;
}
memcpy(lambda, t, (NROOTS + 1) * sizeof(t[0]));
}
}
/* Convert lambda to index form and compute deg(lambda(x)) */
deg_lambda = 0;
for(i = 0; i < NROOTS + 1; i++)
{
lambda[i] = INDEX_OF[lambda[i]];
if(lambda[i] != A0) deg_lambda = i;
}
/* Find roots of the error+erasure locator polynomial by Chien search */
memcpy(&reg[1], &lambda[1], NROOTS * sizeof(reg[0]));
count = 0; /* Number of roots of lambda(x) */
for(i = 1, k = IPRIM - 1; i <= NN; i++, k = MODNN(k + IPRIM))
{
q = 1; /* lambda[0] is always 0 */
for(j = deg_lambda; j > 0; j--)
{
if(reg[j] != A0)
{
reg[j] = MODNN(reg[j] + j);
q ^= ALPHA_TO[reg[j]];
}
}
if (q != 0) continue; /* Not a root */
/* store root (index-form) and error location number */
root[count] = i;
loc[count] = k;
/* If we've already found max possible roots,
* abort the search to save time
*/
if(++count == deg_lambda) break;
}
if(deg_lambda != count)
{
/*
* deg(lambda) unequal to number of roots => uncorrectable
* error detected
*/
count = -1;
goto finish;
}
/*
* Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
* x**NROOTS). in index form. Also find deg(omega).
*/
deg_omega = deg_lambda - 1;
for(i = 0; i <= deg_omega; i++)
{
tmp = 0;
for(j = i; j >= 0; j--)
{
if((s[i - j] != A0) && (lambda[j] != A0))
tmp ^= ALPHA_TO[MODNN(s[i - j] + lambda[j])];
}
omega[i] = INDEX_OF[tmp];
}
/*
* Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
* inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form
*/
for(j = count - 1; j >= 0; j--)
{
num1 = 0;
for(i = deg_omega; i >= 0; i--)
{
if(omega[i] != A0) num1 ^= ALPHA_TO[MODNN(omega[i] + i * root[j])];
}
num2 = ALPHA_TO[MODNN(root[j] * (FCR - 1) + NN)];
den = 0;
/* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
for(i = MIN(deg_lambda, NROOTS - 1) & ~1; i >= 0; i -= 2)
{
if(lambda[i + 1] != A0) den ^= ALPHA_TO[MODNN(lambda[i + 1] + i * root[j])];
}
/* Apply error to data */
if(num1 != 0 && loc[j] >= pad)
{
data[loc[j] - pad] ^= ALPHA_TO[MODNN(INDEX_OF[num1] + INDEX_OF[num2] + NN - INDEX_OF[den])];
}
}
finish:
if(eras_pos != NULL)
{
for(i = 0; i < count; i++) eras_pos[i] = loc[i];
}
return(count);
}