C#でのReed Solomon(リード ソロモン符号)の実装のコードの紹介です。
概要
Reed Solomon(リード ソロモン符号)を実装します。
今回は実装のみで理論の解説はありません。
また、コードも流用コードを含みます。
Delphi版の実装を移植したものとなります
使い方は
こちら
参考・出展
こちらを参考にしております。理論などの解説も下記書籍からどうぞ。
コード
こちらのデモを動作させる場合は、下記 iPentecReedSolomon.cs のみで実行できます。
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace iPentec.ReedSolomonCodec
{
public class iPentecReedSolomon
{
//public static object[] pPP = new object[16 + 1];
public static byte[][] pPP = new byte[16 + 1][];
public const int DEFAULT_MM = 8;
public const int DEFAULT_NP = 32;
// ガロア体原始多項式
public static byte[] PP2 = { 1, 1, 1 };
public static byte[] PP3 = { 1, 1, 0, 1 }; // 1 + x + x^3
public static byte[] PP4 = { 1, 1, 0, 0, 1 }; // 1 + x + x^4
public static byte[] PP5 = { 1, 0, 1, 0, 0, 1 }; // 1 + x^2 + x^5
public static byte[] PP6 = { 1, 1, 0, 0, 0, 0, 1 }; // 1 + x + x^6
public static byte[] PP7 = { 1, 0, 0, 1, 0, 0, 0, 1 }; // 1 + x^3 + x^7
public static byte[] PP8 = { 1, 0, 1, 1, 1, 0, 0, 0, 1 }; // 1+x^2+x^3+x^4+x^8
public static byte[] PP9 = { 1, 0, 0, 0, 1, 0, 0, 0, 0, 1 }; // 1+x^4+x^9
public static byte[] PP10 = { 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1 }; // 1+x^3+x^10
public static byte[] PP11 = { 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x^2+x^11
public static byte[] PP12 = { 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1 }; // 1+x+x^4+x^6+x^12
public static byte[] PP13 = { 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x+x^3+x^4+x^13
public static byte[] PP14 = { 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 }; // 1+x+x^6+x^10+x^14
public static byte[] PP15 = { 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x+x^15
public static byte[] PP16 = { 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1 }; // 1+x+x^3+x^12+x^16
private int FGFDimension = 0;
private int FRSLength = 0;
private int FDataLength = 0;
private int FParityLength = 0;
private int FMaxRSLength = 0;
private byte[] PP = new byte[16 + 1];
// specify irreducible polynomial coeffts
// ガロア体原始多項式
private short[] ExpToVector;
// alpha to
private short[] VectorToExp;
// index of
private short[] gg;
// 誤り訂正生成多項式
private short[] bb;
private short[] data;
private short[] recd;
private byte[] CodeWord;
public int DataLength
{
get
{
return FDataLength;
}
}
public iPentecReedSolomon()
{
pPP[2] = PP2;
pPP[3] = PP3;
pPP[4] = PP4;
pPP[5] = PP5;
pPP[6] = PP6;
pPP[7] = PP7;
pPP[8] = PP8;
pPP[9] = PP9;
pPP[10] = PP10;
pPP[11] = PP11;
pPP[12] = PP12;
pPP[13] = PP13;
pPP[14] = PP14;
pPP[15] = PP15;
pPP[16] = PP16;
InitBuffers();
// Clear the internal buffers
// generate the Galois Field GF(2**mm)
// Generate_gf(DEFAULT_MM);
// FGFDimension:=DEFAULT_MM;
InitParameters(DEFAULT_MM, DEFAULT_NP);
}
public iPentecReedSolomon(int GFDimension, int ParityLength)
{
pPP[2] = PP2;
pPP[3] = PP3;
pPP[4] = PP4;
pPP[5] = PP5;
pPP[6] = PP6;
pPP[7] = PP7;
pPP[8] = PP8;
pPP[9] = PP9;
pPP[10] = PP10;
pPP[11] = PP11;
pPP[12] = PP12;
pPP[13] = PP13;
pPP[14] = PP14;
pPP[15] = PP15;
pPP[16] = PP16;
InitBuffers();
InitParameters(GFDimension, ParityLength);
}
public void dispose()
{
/*
this.Finalize(gg);
this.Finalize(bb);
this.Finalize(data);
this.Finalize(ExpToVector);
this.Finalize(VectorToExp);
this.Finalize(CodeWord);
this.Finalize(recd);
base.Destroy();
*/
}
private void InitParameters(int GFDimension, int ParityLength)
{
FGFDimension = GFDimension;
FParityLength = ParityLength;
FMaxRSLength = (1 << FGFDimension) - 1; // nn=2**mm -1 Max length of codeword
FRSLength = (1 << FGFDimension) - 1; // nn=2**mm -1 Max length of codeword
FDataLength = FRSLength - FParityLength;// data symbols, kk = nn-2*MaxErrors
gg = new short[FParityLength + 1];
bb = new short[FParityLength - 1 + 1];
data = new short[FDataLength + 1];
ExpToVector = new short[FRSLength + 1];
VectorToExp = new short[FRSLength + 1];
CodeWord = new byte[FRSLength - 1 + 1];
recd = new short[FRSLength - 1 + 1];
Generate_GF(FGFDimension);
gen_poly();
}
private void ReInitParameters(int GFDimension, int ParityLength, int RSLength)
{
FGFDimension = GFDimension;
FParityLength = ParityLength;
FMaxRSLength = (1 << FGFDimension) - 1; // nn=2**mm -1 Max length of codeword
FRSLength = RSLength;
FDataLength = FRSLength - FParityLength; // data symbols, kk = nn-2*MaxErrors
gg = new short[FParityLength + 1];
bb = new short[FParityLength - 1 + 1];
data = new short[FDataLength + 1];
ExpToVector = new short[FMaxRSLength + 1];
VectorToExp = new short[FMaxRSLength + 1];
CodeWord = new byte[FRSLength - 1 + 1];
recd = new short[FRSLength - 1 + 1];
Generate_GF(FGFDimension);
gen_poly();
}
protected void Generate_GF(int mm)
{
int i;
short mask;
SetPrimitive(ref PP, mm);
mask = 1;
ExpToVector[mm] = 0;
for (i = 0; i <= (mm - 1); i++) {
ExpToVector[i] = mask;
VectorToExp[ExpToVector[i]] = (short)i;
if (PP[i] != 0) {
ExpToVector[mm] = (short)(ExpToVector[mm] ^ mask);
}
mask = (short)(mask << 1);
}
VectorToExp[ExpToVector[mm]] = (short)mm;
mask = (short)(mask >> 1);
for (i = (mm + 1); i <= (FRSLength - 1); i++) {
if (ExpToVector[i - 1] >= mask) {
ExpToVector[i] = (short)(ExpToVector[mm] ^ ((ExpToVector[i - 1] ^ mask) << 1));
}
else {
ExpToVector[i] = (short)(ExpToVector[i - 1] << 1);
}
VectorToExp[ExpToVector[i]] = (short)i;
}
VectorToExp[0] = -1;
}
protected void gen_poly()
{
short i;
short j;
gg[0] = 2;
// primitive element alpha = 2 for GF(2**mm)
gg[1] = 1;
// g(x) = (X+alpha) initially
for (i = 2; i <= (FRSLength - FDataLength); i++) {
gg[i] = 1;
for (j = (short)(i - 1); j >= 1; j--) {
if ((gg[j] != 0)) {
gg[j] = (short)(gg[j - 1] ^ ExpToVector[(VectorToExp[gg[j]] + i) % FMaxRSLength]);
}
else {
gg[j] = gg[j - 1];
}
}
// gg[0] can never be zero
gg[0] = ExpToVector[(VectorToExp[gg[0]] + i) % FMaxRSLength];
}
// Convert gg[] to index form for quicker encoding.
for (i = 0; i <= FParityLength; i++) {
gg[i] = VectorToExp[gg[i]];
}
}
protected void SetPrimitive(ref byte[] PP, int nIdx)
{
//@ Unsupported function or procedure: 'Move'
// Move(Units.iPentecReedSolomon.pPP[nIdx], PP, (nIdx + 1));
//Array.Copy(pPP, nIdx, PP, 0, nIdx+1);
//Array.Copy(pPP[nIdx], 0, PP, 0, nIdx+1);
for (int i = 0; i < nIdx + 1; i++) {
PP[i] = pPP[nIdx][i];
}
}
// TReedSolomon.InitBuffers
public void InitBuffers()
{
/*
//@ Unsupported function or procedure: 'FillChar'
FillChar(data, sizeof(data), 0);
//@ Unsupported function or procedure: 'FillChar'
FillChar(recd, sizeof(recd), 0);
//@ Unsupported function or procedure: 'FillChar'
FillChar(CodeWord, sizeof(CodeWord), 0);
*/
if (data != null) {
for (int i = 0; i < data.Length; i++) {
data[i] = 0;
}
}
if (recd != null) {
for (int i = 0; i < recd.Length; i++) {
recd[i] = 0;
}
}
if (CodeWord != null) {
for (int i = 0; i < CodeWord.Length; i++) {
CodeWord[i] = 0;
}
}
}
public void EncodeRS(ref byte[] xData, ref byte[] xEncoded)
{
//int e;
byte[] rswork = new byte[123 - 1 + 1];
// RS符号計算用作業領域
int nI;
int i;
int j;
short feedback;
//byte[] axData;
for (i = 0; i < bb.Length; i++) {
bb[i] = 0;
}
for (nI = 0; nI <= (FDataLength - 1); nI++) {
//data[nI] = axData[nI];
data[nI] = xData[nI];
}
for (i = (FDataLength - 1); i >= 0; i--) {
feedback = VectorToExp[data[i] ^ bb[FParityLength - 1]];
if ((feedback != -1)) {
for (j = (FParityLength - 1); j >= 1; j--) {
if ((gg[j] != -1)) {
bb[j] = (short)(bb[j - 1] ^ ExpToVector[(gg[j] + feedback) % FMaxRSLength]);
}
else {
bb[j] = bb[j - 1];
}
}
bb[0] = ExpToVector[(gg[0] + feedback) % FMaxRSLength];
}
else {
for (j = (FParityLength - 1); j >= 1; j--) {
bb[j] = bb[j - 1];
}
bb[0] = 0;
}
}
// put the transmitted codeword, made up of data
// plus parity, in CodeWord
for (nI = 0; nI <= (FParityLength - 1); nI++) {
recd[nI] = bb[nI];
}
for (nI = 0; nI <= (FDataLength - 1); nI++) {
recd[nI + FParityLength] = data[nI];
}
for (nI = 0; nI <= (FRSLength - 1); nI++) {
CodeWord[nI] = (byte)recd[nI];
}
//@ Unsupported function or procedure: 'Move'
//Move(CodeWord[0], xEncoded, FRSLength);
Array.Copy(CodeWord, 0, xEncoded, 0, FRSLength);
}
public void EncodeRSEx(ref byte[] xData, ref byte[] xEncoded, ref byte[] xParity)
{
int nI;
int i;
int j;
short feedback;
int e;
short[] invertgg;
invertgg = new short[gg.Length];
for (i = 0; i < gg.Length; i++) {
//invertgg[i] = gg[gg.Length - 1 - 1 - i];
invertgg[i] = gg[gg.Length - 1 - i];
}
for (i = 0; i < bb.Length; i++) {
bb[i] = 0;
}
//for (nI = 0; nI <= (FRSLength - 1); nI++) {
for (nI = 0; nI <= (FDataLength - 1); nI++) {
data[nI] = xData[nI];
}
for (i = (FDataLength - 1); i >= 0; i--) {
feedback = VectorToExp[data[i] ^ bb[FParityLength - 1]];
if ((feedback != -1)) {
for (j = (FParityLength - 1); j >= 1; j--) {
if ((gg[j] != -1)) {
bb[j] = (short)(bb[j - 1] ^ ExpToVector[(gg[j] + feedback) % FMaxRSLength]);
}
else {
bb[j] = bb[j - 1];
}
}
bb[0] = ExpToVector[(gg[0] + feedback) % FMaxRSLength];
}
else {
for (j = (FParityLength - 1); j >= 1; j--) {
bb[j] = bb[j - 1];
}
bb[0] = 0;
}
}
// put the transmitted codeword, made up of data
// plus parity, in CodeWord
for (nI = 0; nI <= (FParityLength - 1); nI++) {
recd[nI] = bb[nI];
xParity[nI] = (byte)bb[nI];
}
for (nI = 0; nI <= (FDataLength - 1); nI++) {
recd[nI + FParityLength] = data[nI];
}
for (nI = 0; nI <= (FRSLength - 1); nI++) {
CodeWord[nI] = (byte)recd[nI];
}
//@ Unsupported function or procedure: 'Move'
//Move(CodeWord[0], xEncoded, FRSLength);
Array.Copy(CodeWord, 0, xEncoded, 0, FRSLength);
}
//public int DecodeRS(ref object xData, ref object xDecoded)
public int DecodeRS(ref byte[] xData, ref byte[] xDecoded)
{
int result;
//byte[] axData;
//byte[] axDecoded;
//string cStr;
int nI;
int nJ;
int nK;
int i;
int j;
short u;
short q;
short[][] elp; // Array[0..np+1 , 0..np - 1] of smallint ;
short[] d; // Array[0..np+1] of smallint ;
short[] l; // Array[0..np+1] of smallint ;
short[] u_lu; // Array[0..np+1] of smallint ;
short[] s; // Array[0..np] of smallint ;
short count;
short syn_error;
int lpc;
short[] root; // Array[0..MaxErrors - 1] of smallint ;
short[] loc; // Array[0..MaxErrors - 1] of smallint ;
short[] z; // Array[0..MaxErrors] of smallint ;
short[] err; // Array[0..nn-1] of smallint ;
short[] reg; // Array[0..MaxErrors] of smallint ;
// DecodeRS
elp = new short[FParityLength + 1 + 1][];
for (lpc = 0; lpc < elp.Length; lpc++) {
//@ Unsupported function or procedure: 'SetLength'
//elp[lpc].Length = FParityLength - 1 + 1;
elp[lpc] = new short[FParityLength - 1 + 1];
}
d = new short[FParityLength + 1 + 1];
l = new short[FParityLength + 1 + 1];
u_lu = new short[FParityLength + 1 + 1];
s = new short[FParityLength + 1];
root = new short[Convert.ToInt64(FParityLength / 2) - 1 + 1];
loc = new short[Convert.ToInt64(FParityLength / 2) - 1 + 1];
z = new short[Convert.ToInt64(FParityLength / 2) + 1];
err = new short[FRSLength - 1 + 1];
reg = new short[Convert.ToInt64(FParityLength / 2) + 1];
for (nI = 0; nI <= (FRSLength - 1); nI++) {
//recd[nI] = axData[nI];
recd[nI] = xData[nI];
}
for (i = 0; i <= (FRSLength - 1); i++) {
recd[i] = VectorToExp[recd[i]];
}
// put recd[i] into index form
count = 0;
syn_error = 0;
result = 0;
// first form the syndromes
for (i = 1; i <= FParityLength; i++) {
s[i] = 0;
for (j = 0; j <= (FRSLength - 1); j++) {
if (recd[j] != -1) {
// recd[j] in index form
s[i] = (short)(s[i] ^ ExpToVector[(recd[j] + i * j) % FMaxRSLength]);
}
}
// convert syndrome from polynomial form to index form
if ((s[i] != 0)) {
syn_error = 1;
// set flag if non-zero syndrome => error
}
s[i] = VectorToExp[s[i]];
}
if (syn_error != 0) {
d[0] = 0;
// index form
d[1] = s[1];
// index form
elp[0][0] = 0;
// index form
elp[1][0] = 1;
// polynomial form
for (i = 1; i <= (FParityLength - 1); i++) {
elp[0][i] = -1;
// index form
elp[1][i] = 0;
// polynomial form
}
l[0] = 0;
l[1] = 0;
u_lu[0] = -1;
u_lu[1] = 0;
u = 0;
while (((u < FParityLength) && (l[u + 1] <= Convert.ToInt64(FParityLength / 2)))) {
u++;
if ((d[u] == -1)) {
l[u + 1] = l[u];
for (i = 0; i <= l[u]; i++) {
elp[u + 1][i] = elp[u][i];
elp[u][i] = VectorToExp[elp[u][i]];
}
}
else {
// search for words with greatest u_lu[q] for which d[q]!=0
q = (short)(u - 1);
while (((d[q] == -1) && (q > 0))) {
q -= 1;
}
// have found first non-zero d[q]
if ((q > 0)) {
j = q;
while (j > 0) {
j -= 1;
if (((d[j] != -1) && (u_lu[q] < u_lu[j]))) {
q = (short)j;
}
}
}
// have now found q such that d[u]!=0 and u_lu[q] is maximum
// store degree of new elp polynomial
if ((l[u] > l[q] + u - q)) {
l[u + 1] = l[u];
}
else {
l[u + 1] = (short)(l[q] + u - q);
}
// form new elp(x)
for (i = 0; i <= (FParityLength - 1); i++) {
elp[u + 1][i] = 0;
}
for (i = 0; i <= l[q]; i++) {
// if (elp[q][i] <> -1) then elp[u+1][i+u-q] := ExpToVector[(d[u] + FRSLength - d[q] + elp[q][i]) mod FMaxRSLength] ;
if ((elp[q][i] != -1)) {
elp[u + 1][i + u - q] = ExpToVector[(d[u] + FMaxRSLength - d[q] + elp[q][i]) % FMaxRSLength];
}
}
for (i = 0; i <= l[u]; i++) {
elp[u + 1][i] = (short)(elp[u + 1][i] ^ elp[u][i]);
// convert old elp value to index
elp[u][i] = VectorToExp[elp[u][i]];
}
}
u_lu[u + 1] = (short)(u - l[u + 1]);
// form (u+1)th discrepancy
if (u < FParityLength) {
// no discrepancy computed on last iteration
if ((s[u + 1] != -1)) {
d[u + 1] = ExpToVector[s[u + 1]];
}
else {
d[u + 1] = 0;
}
for (i = 1; i <= l[u + 1]; i++) {
if (((s[u + 1 - i] != -1) && (elp[u + 1][i] != 0))) {
d[u + 1] = (short)(d[u + 1] ^ ExpToVector[(s[u + 1 - i] + VectorToExp[elp[u + 1][i]]) % FMaxRSLength]);
}
}
// put d[u+1] into index form
d[u + 1] = VectorToExp[d[u + 1]];
}
}
// end While
u++;
if (l[u] <= Convert.ToInt64(FParityLength / 2)) {
// can correct error
// put elp into index form
for (i = 0; i <= l[u]; i++) {
elp[u][i] = VectorToExp[elp[u][i]];
}
// find roots of the error location polynomial
for (i = 1; i <= l[u]; i++) {
reg[i] = elp[u][i];
}
// for i := 1 to FRSLength do begin
for (i = 1; i <= FMaxRSLength; i++) {
q = 1;
for (j = 1; j <= l[u]; j++) {
if (reg[j] != -1) {
reg[j] = (short)((reg[j] + j) % FMaxRSLength);
q = (short)(q ^ ExpToVector[reg[j]]);
}
}
if (q == 0) {
// store root and error location number indices
root[count] = (short)i;
// loc[count] := FRSLength - i ;
loc[count] = (short)(FMaxRSLength - i);
count++;
}
}
if (count == l[u]) {
// no. roots = degree of elp hence <= tt errors
result = count;
// form polynomial z(x)
for (i = 1; i <= l[u]; i++) {
// Z[0] = 1 always - do not need
if (((s[i] != -1) && (elp[u][i] != -1))) {
z[i] = (short)(ExpToVector[s[i]] ^ ExpToVector[elp[u][i]]);
}
else if (((s[i] != -1) && (elp[u][i] == -1))) {
z[i] = ExpToVector[s[i]];
}
else if (((s[i] == -1) && (elp[u][i] != -1))) {
z[i] = ExpToVector[elp[u][i]];
}
else {
z[i] = 0;
}
for (j = 1; j <= (i - 1); j++) {
if (((s[j] != -1) && (elp[u][i - j] != -1))) {
z[i] = (short)(z[i] ^ ExpToVector[(elp[u][i - j] + s[j]) % FMaxRSLength]);
}
}
// put into index form
z[i] = VectorToExp[z[i]];
}
// evaluate errors at locations given by
// error location numbers loc[i]
for (i = 0; i <= (FRSLength - 1); i++) {
// for i := 0 to (FMaxRSLength - 1) do begin
err[i] = 0;
// convert recd[] to polynomial form
if (recd[i] != -1) {
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
}
for (i = 0; i <= (l[u] - 1); i++) {
// compute numerator of error term first
err[loc[i]] = 1;
// accounts for z[0]
for (j = 1; j <= l[u]; j++) {
if (z[j] != -1) {
err[loc[i]] = (short)(err[loc[i]] ^ ExpToVector[(z[j] + j * root[i]) % FMaxRSLength]);
}
}
if (err[loc[i]] != 0) {
err[loc[i]] = VectorToExp[err[loc[i]]];
q = 0;
// form denominator of error term
for (j = 0; j <= (l[u] - 1); j++) {
if (j != i) {
q = (short)(q + VectorToExp[1 ^ ExpToVector[(loc[j] + root[i]) % FMaxRSLength]]);
}
}
q = (short)(q % FMaxRSLength);
// err[loc[i]] := ExpToVector[(err[loc[i]] - q + FRSLength) mod FMaxRSLength] ;
err[loc[i]] = ExpToVector[(err[loc[i]] - q + FMaxRSLength) % FMaxRSLength];
// recd[i] must be in polynomial form
recd[loc[i]] = (short)(recd[loc[i]] ^ err[loc[i]]);
}
}
}
else {
// no. roots != degree of elp => >tt errors and cannot solve
result = -1;
// Signal an error.
for (i = 0; i <= (FRSLength - 1); i++) {
// could return error flag if desired
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
// just output received codeword as is
}
}
}
else {
// if l[u] <= tt then
// elp has degree has degree >tt hence cannot solve
for (i = 0; i <= (FRSLength - 1); i++) {
// could return error flag if desired
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
// just output received codeword as is
}
}
}
else {
// If syn_error <> 0 then
// no non-zero syndromes => no errors: output received codeword
for (i = 0; i <= (FRSLength - 1); i++) {
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
result = 0;
// No errors ocurred.
}
}
//xDecoded = new byte[FRSLength - 1];
for (nI = 0; nI <= (FRSLength - 1); nI++) {
xDecoded[nI] = (byte)recd[nI];
}
return result;
}
public int DecodeRSEx(ref byte[] xData, ref byte[] xParity, ref byte[] xDecoded)
{
int result;
string cStr;
int nI;
int nJ;
int nK;
int i;
int j;
short u;
short q;
short[][] elp; // Array[0..np+1 , 0..np - 1] of smallint ;
short[] d; // Array[0..np+1] of smallint ;
short[] l; // Array[0..np+1] of smallint ;
short[] u_lu; // Array[0..np+1] of smallint ;
short[] s; // Array[0..np] of smallint ;
short count;
short syn_error;
int lpc;
short[] root; // Array[0..MaxErrors - 1] of smallint ;
short[] loc; // Array[0..MaxErrors - 1] of smallint ;
short[] z; // Array[0..MaxErrors] of smallint ;
short[] err; // Array[0..nn-1] of smallint ;
short[] reg; // Array[0..MaxErrors] of smallint ;
// DecodeRS
elp = new short[FParityLength + 1 + 1][];
for (lpc = 0; lpc < elp.Length; lpc++) {
elp[lpc] = new short[FParityLength - 1 + 1];
//@ Unsupported function or procedure: 'SetLength'
//elp[lpc].Length = FParityLength - 1 + 1;
}
d = new short[FParityLength + 1 + 1];
l = new short[FParityLength + 1 + 1];
u_lu = new short[FParityLength + 1 + 1];
s = new short[FParityLength + 1];
root = new short[Convert.ToInt64(FParityLength / 2) - 1 + 1];
loc = new short[Convert.ToInt64(FParityLength / 2) - 1 + 1];
z = new short[Convert.ToInt64(FParityLength / 2) + 1];
err = new short[FRSLength - 1 + 1];
reg = new short[Convert.ToInt64(FParityLength / 2) + 1];
for (nI = 0; nI <= (FParityLength - 1); nI++) {
recd[nI] = xParity[nI];
}
//for (nI = 0; nI <= (FRSLength - 1); nI++) {
for (nI = 0; nI <= (FDataLength - 1); nI++) {
recd[FParityLength + nI] = xData[nI];
}
for (i = 0; i <= (FRSLength - 1); i++) {
recd[i] = VectorToExp[recd[i]];
}
// put recd[i] into index form
count = 0;
syn_error = 0;
result = 0;
// first form the syndromes
for (i = 1; i <= FParityLength; i++) {
s[i] = 0;
for (j = 0; j <= (FRSLength - 1); j++) {
if (recd[j] != -1) {
// recd[j] in index form
s[i] = (short)(s[i] ^ ExpToVector[(recd[j] + i * j) % FMaxRSLength]);
}
}
// convert syndrome from polynomial form to index form
if ((s[i] != 0)) {
syn_error = 1;
// set flag if non-zero syndrome => error
}
s[i] = VectorToExp[s[i]];
}
if (syn_error != 0) {
// if errors, try and correct
// Compute the error location polynomial via the Berlekamp
// iterative algorithm, following the terminology of Lin and
// Costello : d[u] is the 'mu'th discrepancy, where u='mu'+1
// and 'mu' (the Greek letter!) is the step number ranging from
// -1 to 2*tt (see L&C), l[u] is the degree of the elp at that
// step, and u_l[u] is the difference between the step number
// and the degree of the elp.
// Initialize table entries
d[0] = 0;
// index form
d[1] = s[1];
// index form
elp[0][0] = 0;
// index form
elp[1][0] = 1;
// polynomial form
for (i = 1; i <= (FParityLength - 1); i++) {
elp[0][i] = -1;
// index form
elp[1][i] = 0;
// polynomial form
}
l[0] = 0;
l[1] = 0;
u_lu[0] = -1;
u_lu[1] = 0;
u = 0;
while (((u < FParityLength) && (l[u + 1] <= Convert.ToInt64(FParityLength / 2)))) {
u++;
if ((d[u] == -1)) {
l[u + 1] = l[u];
for (i = 0; i <= l[u]; i++) {
elp[u + 1][i] = elp[u][i];
elp[u][i] = VectorToExp[elp[u][i]];
}
}
else {
// search for words with greatest u_lu[q] for which d[q]!=0
q = (short)(u - 1);
while (((d[q] == -1) && (q > 0))) {
q -= 1;
}
// have found first non-zero d[q]
if ((q > 0)) {
j = q;
while (j > 0) {
j -= 1;
if (((d[j] != -1) && (u_lu[q] < u_lu[j]))) {
q = (short)j;
}
}
}
// have now found q such that d[u]!=0 and u_lu[q] is maximum
// store degree of new elp polynomial
if ((l[u] > l[q] + u - q)) {
l[u + 1] = l[u];
}
else {
l[u + 1] = (short)(l[q] + u - q);
}
// form new elp(x)
for (i = 0; i <= (FParityLength - 1); i++) {
elp[u + 1][i] = 0;
}
for (i = 0; i <= l[q]; i++) {
// if (elp[q][i] <> -1) then elp[u+1][i+u-q] := ExpToVector[(d[u] + FRSLength - d[q] + elp[q][i]) mod FMaxRSLength] ;
if ((elp[q][i] != -1)) {
elp[u + 1][i + u - q] = ExpToVector[(d[u] + FMaxRSLength - d[q] + elp[q][i]) % FMaxRSLength];
}
}
for (i = 0; i <= l[u]; i++) {
elp[u + 1][i] = (short)(elp[u + 1][i] ^ elp[u][i]);
// convert old elp value to index
elp[u][i] = VectorToExp[elp[u][i]];
}
}
u_lu[u + 1] = (short)(u - l[u + 1]);
// form (u+1)th discrepancy
if (u < FParityLength) {
// no discrepancy computed on last iteration
if ((s[u + 1] != -1)) {
d[u + 1] = ExpToVector[s[u + 1]];
}
else {
d[u + 1] = 0;
}
for (i = 1; i <= l[u + 1]; i++) {
if (((s[u + 1 - i] != -1) && (elp[u + 1][i] != 0))) {
d[u + 1] = (short)(d[u + 1] ^ ExpToVector[(s[u + 1 - i] + VectorToExp[elp[u + 1][i]]) % FMaxRSLength]);
}
}
// put d[u+1] into index form
d[u + 1] = VectorToExp[d[u + 1]];
}
}
// end While
u++;
if (l[u] <= Convert.ToInt64(FParityLength / 2)) {
// can correct error
// put elp into index form
for (i = 0; i <= l[u]; i++) {
elp[u][i] = VectorToExp[elp[u][i]];
}
// find roots of the error location polynomial
for (i = 1; i <= l[u]; i++) {
reg[i] = elp[u][i];
}
// for i := 1 to FRSLength do begin
for (i = 1; i <= FMaxRSLength; i++) {
q = 1;
for (j = 1; j <= l[u]; j++) {
if (reg[j] != -1) {
reg[j] = (short)((reg[j] + j) % FMaxRSLength);
q = (short)(q ^ ExpToVector[reg[j]]);
}
}
if (q == 0) {
// store root and error location number indices
root[count] = (short)i;
// loc[count] := FRSLength - i ;
loc[count] = (short)(FMaxRSLength - i);
count++;
}
}
if (count == l[u]) {
// no. roots = degree of elp hence <= tt errors
result = count;
// form polynomial z(x)
for (i = 1; i <= l[u]; i++) {
// Z[0] = 1 always - do not need
if (((s[i] != -1) && (elp[u][i] != -1))) {
z[i] = (short)(ExpToVector[s[i]] ^ ExpToVector[elp[u][i]]);
}
else if (((s[i] != -1) && (elp[u][i] == -1))) {
z[i] = ExpToVector[s[i]];
}
else if (((s[i] == -1) && (elp[u][i] != -1))) {
z[i] = ExpToVector[elp[u][i]];
}
else {
z[i] = 0;
}
for (j = 1; j <= (i - 1); j++) {
if (((s[j] != -1) && (elp[u][i - j] != -1))) {
z[i] = (short)(z[i] ^ ExpToVector[(elp[u][i - j] + s[j]) % FMaxRSLength]);
}
}
// put into index form
z[i] = VectorToExp[z[i]];
}
// evaluate errors at locations given by
// error location numbers loc[i]
for (i = 0; i <= (FRSLength - 1); i++) {
// for i := 0 to (FMaxRSLength - 1) do begin
err[i] = 0;
// convert recd[] to polynomial form
if (recd[i] != -1) {
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
}
for (i = 0; i <= (l[u] - 1); i++) {
// compute numerator of error term first
err[loc[i]] = 1;
// accounts for z[0]
for (j = 1; j <= l[u]; j++) {
if (z[j] != -1) {
err[loc[i]] = (short)(err[loc[i]] ^ ExpToVector[(z[j] + j * root[i]) % FMaxRSLength]);
}
}
if (err[loc[i]] != 0) {
err[loc[i]] = VectorToExp[err[loc[i]]];
q = 0;
// form denominator of error term
for (j = 0; j <= (l[u] - 1); j++) {
if (j != i) {
q = (short)(q + VectorToExp[1 ^ ExpToVector[(loc[j] + root[i]) % FMaxRSLength]]);
}
}
q = (short)(q % FMaxRSLength);
// err[loc[i]] := ExpToVector[(err[loc[i]] - q + FRSLength) mod FMaxRSLength] ;
err[loc[i]] = ExpToVector[(err[loc[i]] - q + FMaxRSLength) % FMaxRSLength];
// recd[i] must be in polynomial form
recd[loc[i]] = (short)(recd[loc[i]] ^ err[loc[i]]);
}
}
}
else {
// no. roots != degree of elp => >tt errors and cannot solve
result = -1;
// Signal an error.
for (i = 0; i <= (FRSLength - 1); i++) {
// could return error flag if desired
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
// just output received codeword as is
}
}
}
else {
// if l[u] <= tt then
// elp has degree has degree >tt hence cannot solve
for (i = 0; i <= (FRSLength - 1); i++) {
// could return error flag if desired
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
// just output received codeword as is
}
}
}
else {
// If syn_error <> 0 then
// no non-zero syndromes => no errors: output received codeword
for (i = 0; i <= (FRSLength - 1); i++) {
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
result = 0;
// No errors ocurred.
}
}
// for nI := FParityLength to (FRSLength - 1) do xDecoded[nI-FParityLength] := Recd[nI] ;
for (nI = FParityLength; nI <= (FRSLength - 1); nI++) {
xDecoded[nI - FParityLength] = (byte)recd[nI];
}
return result;
}
public void SetGFDimension(int value)
{
FGFDimension = value;
//ガロア体の次元が変わったらすべて作成しなおす必要がある
ReInitParameters(FGFDimension, FParityLength, FRSLength);
}
}
}
別実装のコードです。iPentecGalois.cs と iPentecReedSolomonCodec.cs のコードで動作します。
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace iPentec.ReedSolomonCodec
{
public class iPentecGalois
{
private byte[] PP;
// specify irreducible polynomial coeffts
// ガロア体原始多項式
private byte[][] pPP = new byte[16+1][];
private int FTableLength = 0;
public int[] gexp;
public int[] glog;
// ガロア体原始多項式
public static byte[] PP2 = { 1, 1, 1 };
public static byte[] PP3 = { 1, 1, 0, 1 }; // 1 + x + x^3
public static byte[] PP4 = { 1, 1, 0, 0, 1 }; // 1 + x + x^4
public static byte[] PP5 = { 1, 0, 1, 0, 0, 1 }; // 1 + x^2 + x^5
public static byte[] PP6 = { 1, 1, 0, 0, 0, 0, 1 }; // 1 + x + x^6
public static byte[] PP7 = { 1, 0, 0, 1, 0, 0, 0, 1 }; // 1 + x^3 + x^7
public static byte[] PP8 = { 1, 0, 1, 1, 1, 0, 0, 0, 1 }; // 1+x^2+x^3+x^4+x^8
public static byte[] PP9 = { 1, 0, 0, 0, 1, 0, 0, 0, 0, 1 }; // 1+x^4+x^9
public static byte[] PP10 = { 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1 }; // 1+x^3+x^10
public static byte[] PP11 = { 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x^2+x^11
public static byte[] PP12 = { 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1 }; // 1+x+x^4+x^6+x^12
public static byte[] PP13 = { 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x+x^3+x^4+x^13
public static byte[] PP14 = { 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 }; // 1+x+x^6+x^10+x^14
public static byte[] PP15 = { 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x+x^15
public static byte[] PP16 = { 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1 };
// gexp: array [0..512] of integer;
// glog: array [0..512] of integer;
// 1+x+x^3+x^12+x^16
//Constructor Create()
//public iPentecGalois : base()
public iPentecGalois()
{
pPP[2] = PP2;
pPP[3] = PP3;
pPP[4] = PP4;
pPP[5] = PP5;
pPP[6] = PP6;
pPP[7] = PP7;
pPP[8] = PP8;
pPP[9] = PP9;
pPP[10] = PP10;
pPP[11] = PP11;
pPP[12] = PP12;
pPP[13] = PP13;
pPP[14] = PP14;
pPP[15] = PP15;
pPP[16] = PP16;
}
public void InitTables(int Dimension)
{
// * initialize the table of powers of alpha */
// init_exp_table();
init_exp_table(Dimension);
}
private void init_exp_table(int Dimension)
{
int i;
short mask;
FTableLength = (1 << Dimension);
gexp = new int[FTableLength];
glog = new int[FTableLength];
SetPrimitive(ref PP, Dimension);
mask = 1;
gexp[Dimension] = 0;
for (i = 0; i <= (Dimension - 1); i++) {
gexp[i] = mask;
glog[gexp[i]] = i;
if (PP[i] != 0) {
gexp[Dimension] = gexp[Dimension] ^ mask;
}
mask = (short)(mask << 1);
}
glog[gexp[Dimension]] = Dimension;
mask = (short)(mask >> 1);
for (i = (Dimension + 1); i <= (FTableLength - 1); i++) {
if (gexp[i - 1] >= mask) {
gexp[i] = gexp[Dimension] ^ ((gexp[i - 1] ^ mask) << 1);
}
else {
gexp[i] = gexp[i - 1] << 1;
}
glog[gexp[i]] = i;
}
}
private void SetPrimitive(ref byte[] PP, int nIdx)
{
Array.Copy(PP, nIdx, PP, 0, (nIdx + 1));
//@ Unsupported function or procedure: 'Move'
//Move(PP[nIdx], PP, (nIdx + 1));
}
// /* multiplication using logarithms */
public int gmult(int a, int b)
{
int result;
int i;
int j;
if ((a == 0) || (b == 0)) {
result = 0;
}
else {
i = glog[a];
j = glog[b];
// result:=gexp[i+j]; //テーブルの長さが２倍のとき
result = gexp[(i + j) % (FTableLength - 1)];
// テーブルの長さがぴったりのとき
}
return result;
}
public int ginv(int elt)
{
int result;
result = gexp[255 - glog[elt]];
return result;
}
}
}
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace iPentec.ReedSolomonCodec
{
public class iPentecReedSolomonCodec
{
private int NPAR = 0;
private int MAXDEG = 0;
private int FDimension = 0;
private iPentecGalois GT = null;
private int FNErasures = 0;
private int[] ErasureLocs;
private int[] ErrorLocs;
private int NErrors = 0;
private int[] Lambda;
// [MAXDEG];
private int[] Omega;
public int[] pBytes;
public int[] synBytes;
public int[] genPoly;
public int Dimension
{
get
{
return FDimension;
}
set
{
// [MAXDEG];
if (FDimension != value) {
FDimension = value;
}
}
}
//public TiPentecReedSolomonCodec(Component Aowner) : base(Aowner)
public iPentecReedSolomonCodec()
{
FDimension = 8;
GT = new iPentecGalois();
}
/*
protected override void Dispose()
{
}
*/
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// 初期化
public void initialize_ecc(int ParityLength)
{
NPAR = ParityLength;
MAXDEG = NPAR * 2;
// * Initialize the galois field arithmetic tables */
GT.InitTables(FDimension);
// init_galois_tables();
// * Compute the encoder generator polynomial */
genPoly = new int[MAXDEG * 2];
compute_genpoly(NPAR, out genPoly);
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// Decode
//デコード関数
//デコード関数を呼び出したのち check_syndrome でシンドロームを求める
public void decode_data(ref char[] data, int nbytes)
{
int i;
int j;
int sum;
synBytes = new int[MAXDEG];
for (j = 0; j < NPAR; j++) {
// for j:=0 to 8-1 do begin
sum = 0;
for (i = 0; i < nbytes; i++) {
sum = ((byte)data[i]) ^ GT.gmult(GT.gexp[j + 1], sum);
}
synBytes[j] = sum;
}
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// Check Syndrome
// シンドローム算出
// Decodeを事前に呼び出す必要がある
public int check_syndrome()
{
int result;
int i;
int nz;
nz = 0;
for (i = 0; i < NPAR; i++) {
if (synBytes[i] != 0) {
nz = 1;
}
}
result = nz;
return result;
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// Encode
// エンコード関数 コードワード版
public void encode_data(char[] msg, int nbytes, out char[] dst)
{
int i;
int dbyte;
int j;
int[] LFSR;
LFSR = new int[NPAR + 1];
for (i = 0; i <= NPAR; i++) {
LFSR[i] = 0;
}
for (i = 0; i < nbytes; i++) {
dbyte = ((byte)msg[i]) ^ LFSR[NPAR - 1];
for (j = NPAR - 1; j >= 1; j--) {
LFSR[j] = LFSR[j - 1] ^ GT.gmult(genPoly[j], dbyte);
}
LFSR[0] = GT.gmult(genPoly[0], dbyte);
}
pBytes = new int[NPAR];
for (i = 0; i < NPAR; i++) {
pBytes[i] = LFSR[i];
}
build_codeword(msg, nbytes, out dst);
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// Encode Parity
// エンコード関数 パリティ分離版
public void encode_data_parity(char[] msg, int nbytes, out char[] Parity)
{
int i;
int dbyte;
int j;
int[] LFSR;
LFSR = new int[NPAR + 1];
for (i = 0; i <= NPAR; i++) {
LFSR[i] = 0;
}
for (i = 0; i < nbytes; i++) {
dbyte = ((byte)msg[i]) ^ LFSR[NPAR - 1];
for (j = NPAR - 1; j >= 1; j--) {
LFSR[j] = LFSR[j - 1] ^ GT.gmult(genPoly[j], dbyte);
}
LFSR[0] = GT.gmult(genPoly[0], dbyte);
}
Parity = new char[NPAR];
pBytes = new int[NPAR];
for (i = 0; i < NPAR; i++) {
Parity[i] = ((char)LFSR[i]);
}
// for i := 0 to NPAR-1 do Pareity := pBytes[i];
}
public short crc_ccitt(string msg, int leng)
{
short result = 0;
return result;
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// correct errors erasures
// 誤り訂正処理
public int correct_errors_erasures(ref char[] codeword, int csize, int nerasures, int[] erasures)
{
int result;
int r;
int i;
int j;
byte err;
int num;
int denom;
FNErasures = nerasures;
for (i = 0; i < nerasures; i++) {
ErasureLocs[i] = erasures[i];
}
Modified_Berlekamp_Massey();
Find_Roots();
if ((NErrors <= NPAR) && (NErrors > 0)) {
// /* first check for illegal error locs */
for (r = 0; r < NErrors; r++) {
if (ErrorLocs[r] >= csize) {
result = 0;
return result;
}
}
for (r = 0; r < NErrors; r++) {
i = ErrorLocs[r];
// /* evaluate Omega at alpha^(-i) */
num = 0;
for (j = 0; j < MAXDEG; j++) {
num = num ^ GT.gmult(Omega[j], GT.gexp[((255 - i) * j) % 255]);
}
// /* evaluate Lambda' (derivative) at alpha^(-i) ; all odd powers disappear */
denom = 0;
j = 1;
while ((j < MAXDEG)) {
denom = denom ^ GT.gmult(Lambda[j], GT.gexp[((255 - i) * (j - 1)) % 255]);
j += 2;
}
err = (byte)GT.gmult(num, GT.ginv(denom));
//codeword[csize - i - 1] = ((char)((byte)codeword[csize - i - 1]) ^ err);
byte b = (byte)codeword[csize - i - 1];
b = (byte)(b ^ err);
codeword[csize - i - 1] = (char)b;
}
result = 1;
return result;
}
else {
result = 0;
}
return result;
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// Compute GenPoly
// 誤り訂正関数作成
// 誤り訂正関数はエンコード時のみ使用される
protected void compute_genpoly(int nbytes, out int[] genpoly)
{
genpoly = new int[0];
int i;
int[] tp = new int[255 + 1];
int[] tp1 = new int[255 + 1];
// /* multiply (x + a^n) for n = 1 to nbytes */
zero_poly(out tp1);
tp1[0] = 1;
for (i = 1; i <= nbytes; i++) {
zero_poly(out tp);
tp[0] = GT.gexp[i];
// * set up x+a^n */
tp[1] = 1;
mult_polys(out genpoly, tp, tp1);
copy_poly(out tp1, genpoly);
}
// QRコード用 誤り訂正多項式
// for i := 0 to nbytes do begin
// genpoly[i]:=GT.gexp[GT.glog[genpoly[i]] - (nbytes-i)];
// end;
}
protected void build_codeword(char[] msg, int nbytes, out char[] dst)
{
int i;
dst = new char[nbytes];
for (i = 0; i < nbytes; i++) {
dst[i] = msg[i];
}
for (i = 0; i < NPAR; i++) {
dst[i + nbytes] = ((char)pBytes[NPAR - 1 - i]);
}
}
protected void Modified_Berlekamp_Massey()
{
int n;
int l;
int l2;
int k;
int d;
int i;
int[] psi;
int[] psi2;
int[] DA;
int[] gamma;
psi = new int[MAXDEG];
psi2 = new int[MAXDEG];
DA = new int[MAXDEG];
gamma = new int[MAXDEG];
// /* initialize Gamma, the erasure locator polynomial */
init_gamma(out gamma);
// /* initialize to z */
copy_poly(out DA, gamma);
mul_z_poly(ref DA);
copy_poly(out psi, gamma);
k = -1;
l = FNErasures;
for (n = FNErasures; n < NPAR; n++) {
// for n:=FNErasures to 8-1 do begin
d = compute_discrepancy(psi, synBytes, l, n);
if (d != 0) {
// /* psi2 = psi - d*D */
for (i = 0; i < MAXDEG; i++) {
psi2[i] = psi[i] ^ GT.gmult(d, DA[i]);
}
if (l < (n - k)) {
l2 = n - k;
k = n - l;
// /* D = scale_poly(ginv(d), psi); */
for (i = 0; i < MAXDEG; i++) {
DA[i] = GT.gmult(psi[i], GT.ginv(d));
}
l = l2;
}
// /* psi = psi2 */
for (i = 0; i < MAXDEG; i++) {
psi[i] = psi2[i];
}
}
mul_z_poly(ref DA);
}
// λの初期化
Lambda = new int[MAXDEG];
for (i = 0; i < MAXDEG; i++) {
Lambda[i] = psi[i];
}
compute_modified_omega();
}
// * gamma = product (1-z*a^Ij) for erasure locs Ij */
protected void init_gamma(out int[] gamma)
{
int e;
int[] tmp;
tmp = new int[MAXDEG];
zero_poly(out gamma);
zero_poly(out tmp);
gamma[0] = 1;
for (e = 0; e < FNErasures; e++) {
copy_poly(out tmp, gamma);
scale_poly(GT.gexp[ErasureLocs[e]], out tmp);
mul_z_poly(ref tmp);
add_polys(out gamma, tmp);
}
}
protected int compute_discrepancy(int[] lambda, int[] S, int L, int n)
{
int result;
int i;
int sum;
sum = 0;
for (i = 0; i <= L; i++) {
sum = sum ^ GT.gmult(lambda[i], S[n - i]);
}
result = sum;
return result;
}
protected void compute_modified_omega()
{
int i;
int[] product;
product = new int[MAXDEG * 2];
mult_polys(out product, Lambda, synBytes);
// ωの初期化 (0でクリア後 productで初期値を代入)
Omega = new int[MAXDEG];
zero_poly(out Omega);
for (i = 0; i < NPAR; i++) {
Omega[i] = product[i];
}
//this.Finalize(product);
}
protected void Find_Roots()
{
int sum;
int r;
int k;
NErrors = 0;
for (r = 1; r < 256; r++) {
sum = 0;
// /* evaluate lambda at r */
for (k = 0; k < NPAR + 1; k++) {
sum = sum ^ GT.gmult(GT.gexp[(k * r) % 255], Lambda[k]);
}
if ((sum == 0)) {
ErrorLocs[NErrors] = 255 - r;
NErrors++;
}
}
}
public void zero_poly(out int[] poly)
{
int i;
poly = new int[MAXDEG];
for (i = 0; i < MAXDEG; i++) {
poly[i] = 0;
}
}
public void copy_poly(out int[] dst, int[] src)
{
int i;
dst = new int[MAXDEG];
for (i = 0; i < MAXDEG; i++) {
dst[i] = src[i];
}
}
public void add_polys(out int[] dst, int[] src)
{
int i;
dst = new int[MAXDEG];
for (i = 0; i < MAXDEG; i++) {
dst[i] = dst[i] ^ src[i];
}
}
// /* multiply by z, i.e., shift right by 1 */
protected void mul_z_poly(ref int[] src)
{
int i;
src = new int[MAXDEG - 1];
for (i = MAXDEG - 1; i >= 1; i--) {
src[i] = src[i - 1];
}
src[0] = 0;
}
public void scale_poly(int k, out int[] poly)
{
int i;
poly = new int[MAXDEG];
for (i = 0; i < MAXDEG; i++) {
poly[i] = GT.gmult(k, poly[i]);
}
}
public void mult_polys(out int[] dst, int[] p1, int[] p2)
{
int i;
int j;
int[] tmp1;
tmp1 = new int[MAXDEG * 2];
dst = new int[MAXDEG * 2];
for (i = 0; i < (MAXDEG * 2); i++) {
dst[i] = 0;
}
for (i = 0; i < MAXDEG; i++) {
for (j = MAXDEG; j < (MAXDEG * 2); j++) {
tmp1[j] = 0;
}
// /* scale tmp1 by p1[i] */
for (j = 0; j < MAXDEG; j++) {
tmp1[j] = GT.gmult(p2[j], p1[i]);
}
// /* and mult (shift) tmp1 right by i */
j = (MAXDEG * 2) - 1;
while ((j >= i)) {
tmp1[j] = tmp1[j - i];
j -= 1;
}
for (j = 0; j < i; j++) {
tmp1[j] = 0;
}
// /* add into partial product */
for (j = 0; j < (MAXDEG * 2); j++) {
dst[j] = dst[j] ^ tmp1[j];
}
}
}
} // end TiPentecReedSolomonCodec
}
前のバージョンのコード
iPentecGalois.cs
ガロア体原始多項式実装クラス
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace iPentec.ReedSolomonCodec
{
public class iPentecGalois
{
private byte[] PP;
// specify irreducible polynomial coeffts
// ガロア体原始多項式
private byte[][] pPP = new byte[16+1][];
private int FTableLength = 0;
public int[] gexp;
public int[] glog;
// ガロア体原始多項式
public static byte[] PP2 = { 1, 1, 1 };
public static byte[] PP3 = { 1, 1, 0, 1 }; // 1 + x + x^3
public static byte[] PP4 = { 1, 1, 0, 0, 1 }; // 1 + x + x^4
public static byte[] PP5 = { 1, 0, 1, 0, 0, 1 }; // 1 + x^2 + x^5
public static byte[] PP6 = { 1, 1, 0, 0, 0, 0, 1 }; // 1 + x + x^6
public static byte[] PP7 = { 1, 0, 0, 1, 0, 0, 0, 1 }; // 1 + x^3 + x^7
public static byte[] PP8 = { 1, 0, 1, 1, 1, 0, 0, 0, 1 }; // 1+x^2+x^3+x^4+x^8
public static byte[] PP9 = { 1, 0, 0, 0, 1, 0, 0, 0, 0, 1 }; // 1+x^4+x^9
public static byte[] PP10 = { 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1 }; // 1+x^3+x^10
public static byte[] PP11 = { 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x^2+x^11
public static byte[] PP12 = { 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1 }; // 1+x+x^4+x^6+x^12
public static byte[] PP13 = { 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x+x^3+x^4+x^13
public static byte[] PP14 = { 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 }; // 1+x+x^6+x^10+x^14
public static byte[] PP15 = { 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x+x^15
public static byte[] PP16 = { 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1 };
public iPentecGalois()
{
pPP[2] = PP2;
pPP[3] = PP3;
pPP[4] = PP4;
pPP[5] = PP5;
pPP[6] = PP6;
pPP[7] = PP7;
pPP[8] = PP8;
pPP[9] = PP9;
pPP[10] = PP10;
pPP[11] = PP11;
pPP[12] = PP12;
pPP[13] = PP13;
pPP[14] = PP14;
pPP[15] = PP15;
pPP[16] = PP16;
}
public void InitTables(int Dimension)
{
// * initialize the table of powers of alpha */
// init_exp_table();
init_exp_table(Dimension);
}
private void init_exp_table(int Dimension)
{
int i;
short mask;
FTableLength = (1 << Dimension);
gexp = new int[FTableLength];
glog = new int[FTableLength];
SetPrimitive(ref PP, Dimension);
mask = 1;
gexp[Dimension] = 0;
for (i = 0; i <= (Dimension - 1); i++) {
gexp[i] = mask;
glog[gexp[i]] = i;
if (PP[i] != 0) {
gexp[Dimension] = gexp[Dimension] ^ mask;
}
mask = (short)(mask << 1);
}
glog[gexp[Dimension]] = Dimension;
mask = (short)(mask >> 1);
for (i = (Dimension + 1); i <= (FTableLength - 1); i++) {
if (gexp[i - 1] >= mask) {
gexp[i] = gexp[Dimension] ^ ((gexp[i - 1] ^ mask) << 1);
}
else {
gexp[i] = gexp[i - 1] << 1;
}
glog[gexp[i]] = i;
}
}
private void SetPrimitive(ref byte[] PP, int nIdx)
{
Array.Copy(PP, nIdx, PP, 0, (nIdx + 1));
//@ Unsupported function or procedure: 'Move'
//Move(PP[nIdx], PP, (nIdx + 1));
}
// /* multiplication using logarithms */
public int gmult(int a, int b)
{
int result;
int i;
int j;
if ((a == 0) || (b == 0)) {
result = 0;
}
else {
i = glog[a];
j = glog[b];
// result:=gexp[i+j]; //テーブルの長さが２倍のとき
result = gexp[(i + j) % (FTableLength - 1)];
// テーブルの長さがぴったりのとき
}
return result;
}
public int ginv(int elt)
{
int result;
result = gexp[255 - glog[elt]];
return result;
}
}
}
iPentecReedSolomon.cs
Reed Solomon符号化実装クラス
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace iPentec.ReedSolomonCodec
{
public class iPentecReedSolomon
{
//public static object[] pPP = new object[16 + 1];
public static byte[][] pPP = new byte[16 + 1][];
public const int DEFAULT_MM = 8;
public const int DEFAULT_NP = 32;
// ガロア体原始多項式
public static byte[] PP2 = { 1, 1, 1 };
public static byte[] PP3 = { 1, 1, 0, 1 }; // 1 + x + x^3
public static byte[] PP4 = { 1, 1, 0, 0, 1 }; // 1 + x + x^4
public static byte[] PP5 = { 1, 0, 1, 0, 0, 1 }; // 1 + x^2 + x^5
public static byte[] PP6 = { 1, 1, 0, 0, 0, 0, 1 }; // 1 + x + x^6
public static byte[] PP7 = { 1, 0, 0, 1, 0, 0, 0, 1 }; // 1 + x^3 + x^7
public static byte[] PP8 = { 1, 0, 1, 1, 1, 0, 0, 0, 1 }; // 1+x^2+x^3+x^4+x^8
public static byte[] PP9 = { 1, 0, 0, 0, 1, 0, 0, 0, 0, 1 }; // 1+x^4+x^9
public static byte[] PP10 = { 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1 }; // 1+x^3+x^10
public static byte[] PP11 = { 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x^2+x^11
public static byte[] PP12 = { 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1 }; // 1+x+x^4+x^6+x^12
public static byte[] PP13 = { 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x+x^3+x^4+x^13
public static byte[] PP14 = { 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 }; // 1+x+x^6+x^10+x^14
public static byte[] PP15 = { 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }; // 1+x+x^15
public static byte[] PP16 = { 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1 }; // 1+x+x^3+x^12+x^16
private int FGFDimension = 0;
private int FRSLength = 0;
private int FDataLength = 0;
private int FParityLength = 0;
private int FMaxRSLength = 0;
private byte[] PP = new byte[16 + 1];
// specify irreducible polynomial coeffts
// ガロア体原始多項式
private short[] ExpToVector;
// alpha to
private short[] VectorToExp;
// index of
private short[] gg;
// 誤り訂正生成多項式
private short[] bb;
private short[] data;
private short[] recd;
private byte[] CodeWord;
public int DataLength
{
get
{
return FDataLength;
}
}
public iPentecReedSolomon()
{
pPP[2] = PP2;
pPP[3] = PP3;
pPP[4] = PP4;
pPP[5] = PP5;
pPP[6] = PP6;
pPP[7] = PP7;
pPP[8] = PP8;
pPP[9] = PP9;
pPP[10] = PP10;
pPP[11] = PP11;
pPP[12] = PP12;
pPP[13] = PP13;
pPP[14] = PP14;
pPP[15] = PP15;
pPP[16] = PP16;
InitBuffers();
// Clear the internal buffers
// generate the Galois Field GF(2**mm)
// Generate_gf(DEFAULT_MM);
// FGFDimension:=DEFAULT_MM;
InitParameters(DEFAULT_MM, DEFAULT_NP);
}
public iPentecReedSolomon(int GFDimension, int RSLength, int ParityLength)
{
pPP[2] = PP2;
pPP[3] = PP3;
pPP[4] = PP4;
pPP[5] = PP5;
pPP[6] = PP6;
pPP[7] = PP7;
pPP[8] = PP8;
pPP[9] = PP9;
pPP[10] = PP10;
pPP[11] = PP11;
pPP[12] = PP12;
pPP[13] = PP13;
pPP[14] = PP14;
pPP[15] = PP15;
pPP[16] = PP16;
InitBuffers();
InitParameters(DEFAULT_MM, DEFAULT_NP);
SetGFDimension(GFDimension);
SetRSLength(RSLength, ParityLength);
}
public void dispose()
{
}
private void InitParameters(int GFDimension, int ParityLength)
{
FGFDimension = GFDimension;
FParityLength = ParityLength;
FMaxRSLength = (1 << FGFDimension) - 1; // nn=2**mm -1 Max length of codeword
FRSLength = (1 << FGFDimension) - 1; // nn=2**mm -1 Max length of codeword
FDataLength = FRSLength - FParityLength;// data symbols, kk = nn-2*MaxErrors
gg = new short[FParityLength + 1];
bb = new short[FParityLength - 1 + 1];
data = new short[FDataLength + 1];
ExpToVector = new short[FRSLength + 1];
VectorToExp = new short[FRSLength + 1];
CodeWord = new byte[FRSLength - 1 + 1];
recd = new short[FRSLength - 1 + 1];
Generate_GF(FGFDimension);
gen_poly();
}
private void ReInitParameters(int GFDimension, int ParityLength, int RSLength)
{
FGFDimension = GFDimension;
FParityLength = ParityLength;
FMaxRSLength = (1 << FGFDimension) - 1; // nn=2**mm -1 Max length of codeword
FRSLength = RSLength;
FDataLength = FRSLength - FParityLength; // data symbols, kk = nn-2*MaxErrors
gg = new short[FParityLength + 1];
bb = new short[FParityLength - 1 + 1];
data = new short[FDataLength + 1];
ExpToVector = new short[FMaxRSLength + 1];
VectorToExp = new short[FMaxRSLength + 1];
CodeWord = new byte[FRSLength - 1 + 1];
recd = new short[FRSLength - 1 + 1];
Generate_GF(FGFDimension);
gen_poly();
}
protected void Generate_GF(int mm)
{
int i;
short mask;
SetPrimitive(ref PP, mm);
mask = 1;
ExpToVector[mm] = 0;
for (i = 0; i <= (mm - 1); i++) {
ExpToVector[i] = mask;
VectorToExp[ExpToVector[i]] = (short)i;
if (PP[i] != 0) {
ExpToVector[mm] = (short)(ExpToVector[mm] ^ mask);
}
mask = (short)(mask << 1);
}
VectorToExp[ExpToVector[mm]] = (short)mm;
mask = (short)(mask >> 1);
for (i = (mm + 1); i <= (FRSLength - 1); i++) {
if (ExpToVector[i - 1] >= mask) {
ExpToVector[i] = (short)(ExpToVector[mm] ^ ((ExpToVector[i - 1] ^ mask) << 1));
}
else {
ExpToVector[i] = (short)(ExpToVector[i - 1] << 1);
}
VectorToExp[ExpToVector[i]] = (short)i;
}
VectorToExp[0] = -1;
}
protected void gen_poly()
{
short i;
short j;
gg[0] = 2;
// primitive element alpha = 2 for GF(2**mm)
gg[1] = 1;
// g(x) = (X+alpha) initially
for (i = 2; i <= (FRSLength - FDataLength); i++) {
gg[i] = 1;
for (j = (short)(i - 1); j >= 1; j--) {
if ((gg[j] != 0)) {
gg[j] = (short)(gg[j - 1] ^ ExpToVector[(VectorToExp[gg[j]] + i) % FMaxRSLength]);
}
else {
gg[j] = gg[j - 1];
}
}
// gg[0] can never be zero
gg[0] = ExpToVector[(VectorToExp[gg[0]] + i) % FMaxRSLength];
}
// Convert gg[] to index form for quicker encoding.
for (i = 0; i <= FParityLength; i++) {
gg[i] = VectorToExp[gg[i]];
}
}
protected void SetPrimitive(ref byte[] PP, int nIdx)
{
//@ Unsupported function or procedure: 'Move'
// Move(Units.iPentecReedSolomon.pPP[nIdx], PP, (nIdx + 1));
//Array.Copy(pPP, nIdx, PP, 0, nIdx+1);
//Array.Copy(pPP[nIdx], 0, PP, 0, nIdx+1);
for (int i = 0; i < nIdx + 1; i++) {
PP[i] = pPP[nIdx][i];
}
}
// TReedSolomon.InitBuffers
public void InitBuffers()
{
/*
//@ Unsupported function or procedure: 'FillChar'
FillChar(data, sizeof(data), 0);
//@ Unsupported function or procedure: 'FillChar'
FillChar(recd, sizeof(recd), 0);
//@ Unsupported function or procedure: 'FillChar'
FillChar(CodeWord, sizeof(CodeWord), 0);
*/
if (data != null) {
for (int i = 0; i < data.Length; i++) {
data[i] = 0;
}
}
if (recd != null) {
for (int i = 0; i < recd.Length; i++) {
recd[i] = 0;
}
}
if (CodeWord != null) {
for (int i = 0; i < CodeWord.Length; i++) {
CodeWord[i] = 0;
}
}
}
public void EncodeRS(ref byte[] xData, ref byte[] xEncoded)
{
//int e;
byte[] rswork = new byte[123 - 1 + 1];
// RS符号計算用作業領域
int nI;
int i;
int j;
short feedback;
//byte[] axData;
for (i = 0; i < bb.Length; i++) {
bb[i] = 0;
}
for (nI = 0; nI <= (FDataLength - 1); nI++) {
//data[nI] = axData[nI];
data[nI] = xData[nI];
}
for (i = (FDataLength - 1); i >= 0; i--) {
feedback = VectorToExp[data[i] ^ bb[FParityLength - 1]];
if ((feedback != -1)) {
for (j = (FParityLength - 1); j >= 1; j--) {
if ((gg[j] != -1)) {
bb[j] = (short)(bb[j - 1] ^ ExpToVector[(gg[j] + feedback) % FMaxRSLength]);
}
else {
bb[j] = bb[j - 1];
}
}
bb[0] = ExpToVector[(gg[0] + feedback) % FMaxRSLength];
}
else {
for (j = (FParityLength - 1); j >= 1; j--) {
bb[j] = bb[j - 1];
}
bb[0] = 0;
}
}
// put the transmitted codeword, made up of data
// plus parity, in CodeWord
for (nI = 0; nI <= (FParityLength - 1); nI++) {
recd[nI] = bb[nI];
}
for (nI = 0; nI <= (FDataLength - 1); nI++) {
recd[nI + FParityLength] = data[nI];
}
for (nI = 0; nI <= (FRSLength - 1); nI++) {
CodeWord[nI] = (byte)recd[nI];
}
//@ Unsupported function or procedure: 'Move'
//Move(CodeWord[0], xEncoded, FRSLength);
Array.Copy(CodeWord, 0, xEncoded, 0, FRSLength);
}
public void EncodeRSEx(ref byte[] xData, ref byte[] xEncoded, ref byte[] xParity)
{
int nI;
int i;
int j;
short feedback;
int e;
short[] invertgg;
invertgg = new short[gg.Length];
for (i = 0; i < gg.Length; i++) {
//invertgg[i] = gg[gg.Length - 1 - 1 - i];
invertgg[i] = gg[gg.Length - 1 - i];
}
for (i = 0; i < bb.Length; i++) {
bb[i] = 0;
}
//for (nI = 0; nI <= (FRSLength - 1); nI++) {
for (nI = 0; nI <= (FDataLength - 1); nI++) {
data[nI] = xData[nI];
}
for (i = (FDataLength - 1); i >= 0; i--) {
feedback = VectorToExp[data[i] ^ bb[FParityLength - 1]];
if ((feedback != -1)) {
for (j = (FParityLength - 1); j >= 1; j--) {
if ((gg[j] != -1)) {
bb[j] = (short)(bb[j - 1] ^ ExpToVector[(gg[j] + feedback) % FMaxRSLength]);
}
else {
bb[j] = bb[j - 1];
}
}
bb[0] = ExpToVector[(gg[0] + feedback) % FMaxRSLength];
}
else {
for (j = (FParityLength - 1); j >= 1; j--) {
bb[j] = bb[j - 1];
}
bb[0] = 0;
}
}
// put the transmitted codeword, made up of data
// plus parity, in CodeWord
for (nI = 0; nI <= (FParityLength - 1); nI++) {
recd[nI] = bb[nI];
xParity[nI] = (byte)bb[nI];
}
for (nI = 0; nI <= (FDataLength - 1); nI++) {
recd[nI + FParityLength] = data[nI];
}
for (nI = 0; nI <= (FRSLength - 1); nI++) {
CodeWord[nI] = (byte)recd[nI];
}
//@ Unsupported function or procedure: 'Move'
//Move(CodeWord[0], xEncoded, FRSLength);
Array.Copy(CodeWord, 0, xEncoded, 0, FRSLength);
}
//public int DecodeRS(ref object xData, ref object xDecoded)
public int DecodeRS(ref byte[] xData, ref byte[] xDecoded)
{
int result;
//byte[] axData;
//byte[] axDecoded;
//string cStr;
int nI;
int nJ;
int nK;
int i;
int j;
short u;
short q;
short[][] elp; // Array[0..np+1 , 0..np - 1] of smallint ;
short[] d; // Array[0..np+1] of smallint ;
short[] l; // Array[0..np+1] of smallint ;
short[] u_lu; // Array[0..np+1] of smallint ;
short[] s; // Array[0..np] of smallint ;
short count;
short syn_error;
int lpc;
short[] root; // Array[0..MaxErrors - 1] of smallint ;
short[] loc; // Array[0..MaxErrors - 1] of smallint ;
short[] z; // Array[0..MaxErrors] of smallint ;
short[] err; // Array[0..nn-1] of smallint ;
short[] reg; // Array[0..MaxErrors] of smallint ;
// DecodeRS
elp = new short[FParityLength + 1 + 1][];
for (lpc = 0; lpc < elp.Length; lpc++) {
//@ Unsupported function or procedure: 'SetLength'
//elp[lpc].Length = FParityLength - 1 + 1;
elp[lpc] = new short[FParityLength - 1 + 1];
}
d = new short[FParityLength + 1 + 1];
l = new short[FParityLength + 1 + 1];
u_lu = new short[FParityLength + 1 + 1];
s = new short[FParityLength + 1];
root = new short[Convert.ToInt64(FParityLength / 2) - 1 + 1];
loc = new short[Convert.ToInt64(FParityLength / 2) - 1 + 1];
z = new short[Convert.ToInt64(FParityLength / 2) + 1];
err = new short[FRSLength - 1 + 1];
reg = new short[Convert.ToInt64(FParityLength / 2) + 1];
for (nI = 0; nI <= (FRSLength - 1); nI++) {
//recd[nI] = axData[nI];
recd[nI] = xData[nI];
}
for (i = 0; i <= (FRSLength - 1); i++) {
recd[i] = VectorToExp[recd[i]];
}
// put recd[i] into index form
count = 0;
syn_error = 0;
result = 0;
// first form the syndromes
for (i = 1; i <= FParityLength; i++) {
s[i] = 0;
for (j = 0; j <= (FRSLength - 1); j++) {
if (recd[j] != -1) {
// recd[j] in index form
s[i] = (short)(s[i] ^ ExpToVector[(recd[j] + i * j) % FMaxRSLength]);
}
}
// convert syndrome from polynomial form to index form
if ((s[i] != 0)) {
syn_error = 1;
// set flag if non-zero syndrome => error
}
s[i] = VectorToExp[s[i]];
}
if (syn_error != 0) {
d[0] = 0;
// index form
d[1] = s[1];
// index form
elp[0][0] = 0;
// index form
elp[1][0] = 1;
// polynomial form
for (i = 1; i <= (FParityLength - 1); i++) {
elp[0][i] = -1;
// index form
elp[1][i] = 0;
// polynomial form
}
l[0] = 0;
l[1] = 0;
u_lu[0] = -1;
u_lu[1] = 0;
u = 0;
while (((u < FParityLength) && (l[u + 1] <= Convert.ToInt64(FParityLength / 2)))) {
u++;
if ((d[u] == -1)) {
l[u + 1] = l[u];
for (i = 0; i <= l[u]; i++) {
elp[u + 1][i] = elp[u][i];
elp[u][i] = VectorToExp[elp[u][i]];
}
}
else {
// search for words with greatest u_lu[q] for which d[q]!=0
q = (short)(u - 1);
while (((d[q] == -1) && (q > 0))) {
q -= 1;
}
// have found first non-zero d[q]
if ((q > 0)) {
j = q;
while (j > 0) {
j -= 1;
if (((d[j] != -1) && (u_lu[q] < u_lu[j]))) {
q = (short)j;
}
}
}
// have now found q such that d[u]!=0 and u_lu[q] is maximum
// store degree of new elp polynomial
if ((l[u] > l[q] + u - q)) {
l[u + 1] = l[u];
}
else {
l[u + 1] = (short)(l[q] + u - q);
}
// form new elp(x)
for (i = 0; i <= (FParityLength - 1); i++) {
elp[u + 1][i] = 0;
}
for (i = 0; i <= l[q]; i++) {
// if (elp[q][i] <> -1) then elp[u+1][i+u-q] := ExpToVector[(d[u] + FRSLength - d[q] + elp[q][i]) mod FMaxRSLength] ;
if ((elp[q][i] != -1)) {
elp[u + 1][i + u - q] = ExpToVector[(d[u] + FMaxRSLength - d[q] + elp[q][i]) % FMaxRSLength];
}
}
for (i = 0; i <= l[u]; i++) {
elp[u + 1][i] = (short)(elp[u + 1][i] ^ elp[u][i]);
// convert old elp value to index
elp[u][i] = VectorToExp[elp[u][i]];
}
}
u_lu[u + 1] = (short)(u - l[u + 1]);
// form (u+1)th discrepancy
if (u < FParityLength) {
// no discrepancy computed on last iteration
if ((s[u + 1] != -1)) {
d[u + 1] = ExpToVector[s[u + 1]];
}
else {
d[u + 1] = 0;
}
for (i = 1; i <= l[u + 1]; i++) {
if (((s[u + 1 - i] != -1) && (elp[u + 1][i] != 0))) {
d[u + 1] = (short)(d[u + 1] ^ ExpToVector[(s[u + 1 - i] + VectorToExp[elp[u + 1][i]]) % FMaxRSLength]);
}
}
// put d[u+1] into index form
d[u + 1] = VectorToExp[d[u + 1]];
}
}
// end While
u++;
if (l[u] <= Convert.ToInt64(FParityLength / 2)) {
// can correct error
// put elp into index form
for (i = 0; i <= l[u]; i++) {
elp[u][i] = VectorToExp[elp[u][i]];
}
// find roots of the error location polynomial
for (i = 1; i <= l[u]; i++) {
reg[i] = elp[u][i];
}
// for i := 1 to FRSLength do begin
for (i = 1; i <= FMaxRSLength; i++) {
q = 1;
for (j = 1; j <= l[u]; j++) {
if (reg[j] != -1) {
reg[j] = (short)((reg[j] + j) % FMaxRSLength);
q = (short)(q ^ ExpToVector[reg[j]]);
}
}
if (q == 0) {
// store root and error location number indices
root[count] = (short)i;
// loc[count] := FRSLength - i ;
loc[count] = (short)(FMaxRSLength - i);
count++;
}
}
if (count == l[u]) {
// no. roots = degree of elp hence <= tt errors
result = count;
// form polynomial z(x)
for (i = 1; i <= l[u]; i++) {
// Z[0] = 1 always - do not need
if (((s[i] != -1) && (elp[u][i] != -1))) {
z[i] = (short)(ExpToVector[s[i]] ^ ExpToVector[elp[u][i]]);
}
else if (((s[i] != -1) && (elp[u][i] == -1))) {
z[i] = ExpToVector[s[i]];
}
else if (((s[i] == -1) && (elp[u][i] != -1))) {
z[i] = ExpToVector[elp[u][i]];
}
else {
z[i] = 0;
}
for (j = 1; j <= (i - 1); j++) {
if (((s[j] != -1) && (elp[u][i - j] != -1))) {
z[i] = (short)(z[i] ^ ExpToVector[(elp[u][i - j] + s[j]) % FMaxRSLength]);
}
}
// put into index form
z[i] = VectorToExp[z[i]];
}
// evaluate errors at locations given by
// error location numbers loc[i]
for (i = 0; i <= (FRSLength - 1); i++) {
// for i := 0 to (FMaxRSLength - 1) do begin
err[i] = 0;
// convert recd[] to polynomial form
if (recd[i] != -1) {
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
}
for (i = 0; i <= (l[u] - 1); i++) {
// compute numerator of error term first
err[loc[i]] = 1;
// accounts for z[0]
for (j = 1; j <= l[u]; j++) {
if (z[j] != -1) {
err[loc[i]] = (short)(err[loc[i]] ^ ExpToVector[(z[j] + j * root[i]) % FMaxRSLength]);
}
}
if (err[loc[i]] != 0) {
err[loc[i]] = VectorToExp[err[loc[i]]];
q = 0;
// form denominator of error term
for (j = 0; j <= (l[u] - 1); j++) {
if (j != i) {
q = (short)(q + VectorToExp[1 ^ ExpToVector[(loc[j] + root[i]) % FMaxRSLength]]);
}
}
q = (short)(q % FMaxRSLength);
// err[loc[i]] := ExpToVector[(err[loc[i]] - q + FRSLength) mod FMaxRSLength] ;
err[loc[i]] = ExpToVector[(err[loc[i]] - q + FMaxRSLength) % FMaxRSLength];
// recd[i] must be in polynomial form
recd[loc[i]] = (short)(recd[loc[i]] ^ err[loc[i]]);
}
}
}
else {
// no. roots != degree of elp => >tt errors and cannot solve
result = -1;
// Signal an error.
for (i = 0; i <= (FRSLength - 1); i++) {
// could return error flag if desired
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
// just output received codeword as is
}
}
}
else {
// if l[u] <= tt then
// elp has degree has degree >tt hence cannot solve
for (i = 0; i <= (FRSLength - 1); i++) {
// could return error flag if desired
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
// just output received codeword as is
}
}
}
else {
// If syn_error <> 0 then
// no non-zero syndromes => no errors: output received codeword
for (i = 0; i <= (FRSLength - 1); i++) {
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
result = 0;
// No errors ocurred.
}
}
//xDecoded = new byte[FRSLength - 1];
for (nI = 0; nI <= (FRSLength - 1); nI++) {
xDecoded[nI] = (byte)recd[nI];
}
return result;
}
public int DecodeRSEx(ref byte[] xData, ref byte[] xParity, ref byte[] xDecoded)
{
int result;
string cStr;
int nI;
int nJ;
int nK;
int i;
int j;
short u;
short q;
short[][] elp; // Array[0..np+1 , 0..np - 1] of smallint ;
short[] d; // Array[0..np+1] of smallint ;
short[] l; // Array[0..np+1] of smallint ;
short[] u_lu; // Array[0..np+1] of smallint ;
short[] s; // Array[0..np] of smallint ;
short count;
short syn_error;
int lpc;
short[] root; // Array[0..MaxErrors - 1] of smallint ;
short[] loc; // Array[0..MaxErrors - 1] of smallint ;
short[] z; // Array[0..MaxErrors] of smallint ;
short[] err; // Array[0..nn-1] of smallint ;
short[] reg; // Array[0..MaxErrors] of smallint ;
// DecodeRS
elp = new short[FParityLength + 1 + 1][];
for (lpc = 0; lpc < elp.Length; lpc++) {
elp[lpc] = new short[FParityLength - 1 + 1];
//@ Unsupported function or procedure: 'SetLength'
//elp[lpc].Length = FParityLength - 1 + 1;
}
d = new short[FParityLength + 1 + 1];
l = new short[FParityLength + 1 + 1];
u_lu = new short[FParityLength + 1 + 1];
s = new short[FParityLength + 1];
root = new short[Convert.ToInt64(FParityLength / 2) - 1 + 1];
loc = new short[Convert.ToInt64(FParityLength / 2) - 1 + 1];
z = new short[Convert.ToInt64(FParityLength / 2) + 1];
err = new short[FRSLength - 1 + 1];
reg = new short[Convert.ToInt64(FParityLength / 2) + 1];
for (nI = 0; nI <= (FParityLength - 1); nI++) {
recd[nI] = xParity[nI];
}
//for (nI = 0; nI <= (FRSLength - 1); nI++) {
for (nI = 0; nI <= (FDataLength - 1); nI++) {
recd[FParityLength + nI] = xData[nI];
}
for (i = 0; i <= (FRSLength - 1); i++) {
recd[i] = VectorToExp[recd[i]];
}
// put recd[i] into index form
count = 0;
syn_error = 0;
result = 0;
// first form the syndromes
for (i = 1; i <= FParityLength; i++) {
s[i] = 0;
for (j = 0; j <= (FRSLength - 1); j++) {
if (recd[j] != -1) {
// recd[j] in index form
s[i] = (short)(s[i] ^ ExpToVector[(recd[j] + i * j) % FMaxRSLength]);
}
}
// convert syndrome from polynomial form to index form
if ((s[i] != 0)) {
syn_error = 1;
// set flag if non-zero syndrome => error
}
s[i] = VectorToExp[s[i]];
}
if (syn_error != 0) {
// if errors, try and correct
// Compute the error location polynomial via the Berlekamp
// iterative algorithm, following the terminology of Lin and
// Costello : d[u] is the 'mu'th discrepancy, where u='mu'+1
// and 'mu' (the Greek letter!) is the step number ranging from
// -1 to 2*tt (see L&C), l[u] is the degree of the elp at that
// step, and u_l[u] is the difference between the step number
// and the degree of the elp.
// Initialize table entries
d[0] = 0;
// index form
d[1] = s[1];
// index form
elp[0][0] = 0;
// index form
elp[1][0] = 1;
// polynomial form
for (i = 1; i <= (FParityLength - 1); i++) {
elp[0][i] = -1;
// index form
elp[1][i] = 0;
// polynomial form
}
l[0] = 0;
l[1] = 0;
u_lu[0] = -1;
u_lu[1] = 0;
u = 0;
while (((u < FParityLength) && (l[u + 1] <= Convert.ToInt64(FParityLength / 2)))) {
u++;
if ((d[u] == -1)) {
l[u + 1] = l[u];
for (i = 0; i <= l[u]; i++) {
elp[u + 1][i] = elp[u][i];
elp[u][i] = VectorToExp[elp[u][i]];
}
}
else {
// search for words with greatest u_lu[q] for which d[q]!=0
q = (short)(u - 1);
while (((d[q] == -1) && (q > 0))) {
q -= 1;
}
// have found first non-zero d[q]
if ((q > 0)) {
j = q;
while (j > 0) {
j -= 1;
if (((d[j] != -1) && (u_lu[q] < u_lu[j]))) {
q = (short)j;
}
}
}
// have now found q such that d[u]!=0 and u_lu[q] is maximum
// store degree of new elp polynomial
if ((l[u] > l[q] + u - q)) {
l[u + 1] = l[u];
}
else {
l[u + 1] = (short)(l[q] + u - q);
}
// form new elp(x)
for (i = 0; i <= (FParityLength - 1); i++) {
elp[u + 1][i] = 0;
}
for (i = 0; i <= l[q]; i++) {
// if (elp[q][i] <> -1) then elp[u+1][i+u-q] := ExpToVector[(d[u] + FRSLength - d[q] + elp[q][i]) mod FMaxRSLength] ;
if ((elp[q][i] != -1)) {
elp[u + 1][i + u - q] = ExpToVector[(d[u] + FMaxRSLength - d[q] + elp[q][i]) % FMaxRSLength];
}
}
for (i = 0; i <= l[u]; i++) {
elp[u + 1][i] = (short)(elp[u + 1][i] ^ elp[u][i]);
// convert old elp value to index
elp[u][i] = VectorToExp[elp[u][i]];
}
}
u_lu[u + 1] = (short)(u - l[u + 1]);
// form (u+1)th discrepancy
if (u < FParityLength) {
// no discrepancy computed on last iteration
if ((s[u + 1] != -1)) {
d[u + 1] = ExpToVector[s[u + 1]];
}
else {
d[u + 1] = 0;
}
for (i = 1; i <= l[u + 1]; i++) {
if (((s[u + 1 - i] != -1) && (elp[u + 1][i] != 0))) {
d[u + 1] = (short)(d[u + 1] ^ ExpToVector[(s[u + 1 - i] + VectorToExp[elp[u + 1][i]]) % FMaxRSLength]);
}
}
// put d[u+1] into index form
d[u + 1] = VectorToExp[d[u + 1]];
}
}
// end While
u++;
if (l[u] <= Convert.ToInt64(FParityLength / 2)) {
// can correct error
// put elp into index form
for (i = 0; i <= l[u]; i++) {
elp[u][i] = VectorToExp[elp[u][i]];
}
// find roots of the error location polynomial
for (i = 1; i <= l[u]; i++) {
reg[i] = elp[u][i];
}
// for i := 1 to FRSLength do begin
for (i = 1; i <= FMaxRSLength; i++) {
q = 1;
for (j = 1; j <= l[u]; j++) {
if (reg[j] != -1) {
reg[j] = (short)((reg[j] + j) % FMaxRSLength);
q = (short)(q ^ ExpToVector[reg[j]]);
}
}
if (q == 0) {
// store root and error location number indices
root[count] = (short)i;
// loc[count] := FRSLength - i ;
loc[count] = (short)(FMaxRSLength - i);
count++;
}
}
if (count == l[u]) {
// no. roots = degree of elp hence <= tt errors
result = count;
// form polynomial z(x)
for (i = 1; i <= l[u]; i++) {
// Z[0] = 1 always - do not need
if (((s[i] != -1) && (elp[u][i] != -1))) {
z[i] = (short)(ExpToVector[s[i]] ^ ExpToVector[elp[u][i]]);
}
else if (((s[i] != -1) && (elp[u][i] == -1))) {
z[i] = ExpToVector[s[i]];
}
else if (((s[i] == -1) && (elp[u][i] != -1))) {
z[i] = ExpToVector[elp[u][i]];
}
else {
z[i] = 0;
}
for (j = 1; j <= (i - 1); j++) {
if (((s[j] != -1) && (elp[u][i - j] != -1))) {
z[i] = (short)(z[i] ^ ExpToVector[(elp[u][i - j] + s[j]) % FMaxRSLength]);
}
}
// put into index form
z[i] = VectorToExp[z[i]];
}
// evaluate errors at locations given by
// error location numbers loc[i]
for (i = 0; i <= (FRSLength - 1); i++) {
// for i := 0 to (FMaxRSLength - 1) do begin
err[i] = 0;
// convert recd[] to polynomial form
if (recd[i] != -1) {
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
}
for (i = 0; i <= (l[u] - 1); i++) {
// compute numerator of error term first
err[loc[i]] = 1;
// accounts for z[0]
for (j = 1; j <= l[u]; j++) {
if (z[j] != -1) {
err[loc[i]] = (short)(err[loc[i]] ^ ExpToVector[(z[j] + j * root[i]) % FMaxRSLength]);
}
}
if (err[loc[i]] != 0) {
err[loc[i]] = VectorToExp[err[loc[i]]];
q = 0;
// form denominator of error term
for (j = 0; j <= (l[u] - 1); j++) {
if (j != i) {
q = (short)(q + VectorToExp[1 ^ ExpToVector[(loc[j] + root[i]) % FMaxRSLength]]);
}
}
q = (short)(q % FMaxRSLength);
// err[loc[i]] := ExpToVector[(err[loc[i]] - q + FRSLength) mod FMaxRSLength] ;
err[loc[i]] = ExpToVector[(err[loc[i]] - q + FMaxRSLength) % FMaxRSLength];
// recd[i] must be in polynomial form
recd[loc[i]] = (short)(recd[loc[i]] ^ err[loc[i]]);
}
}
}
else {
// no. roots != degree of elp => >tt errors and cannot solve
result = -1;
// Signal an error.
for (i = 0; i <= (FRSLength - 1); i++) {
// could return error flag if desired
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
// just output received codeword as is
}
}
}
else {
// if l[u] <= tt then
// elp has degree has degree >tt hence cannot solve
for (i = 0; i <= (FRSLength - 1); i++) {
// could return error flag if desired
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
// just output received codeword as is
}
}
}
else {
// If syn_error <> 0 then
// no non-zero syndromes => no errors: output received codeword
for (i = 0; i <= (FRSLength - 1); i++) {
if (recd[i] != -1) {
// convert recd[] to polynomial form
recd[i] = ExpToVector[recd[i]];
}
else {
recd[i] = 0;
}
result = 0;
// No errors ocurred.
}
}
// for nI := FParityLength to (FRSLength - 1) do xDecoded[nI-FParityLength] := Recd[nI] ;
for (nI = FParityLength; nI <= (FRSLength - 1); nI++) {
xDecoded[nI - FParityLength] = (byte)recd[nI];
}
return result;
}
public void SetGFDimension(int value)
{
FGFDimension = value;
//ガロア体の次元が変わったらすべて作成しなおす必要がある
ReInitParameters(FGFDimension, FParityLength, FRSLength);
}
public void SetRSLength(int RSLength, int ParityLength)
{
FRSLength = RSLength;
FParityLength = ParityLength;
FDataLength = FRSLength - FParityLength; //data symbols, kk = nn-2*MaxErrors
gg = new short[FParityLength + 1];
bb = new short[FParityLength];
data = new short[FDataLength + 1];
CodeWord = new byte[FRSLength - 1 + 1];
recd = new short[FRSLength - 1 + 1];
gen_poly();
}
public void SetRSMaxLength()
{
int ml;
ml = (1 << FGFDimension) - 1;
// nn=2**mm -1 length of codeword
//RSLength = ml;
SetRSLength(ml, FParityLength);
}
}
} // end TiPentecReedSolomon
iPentecReedSolomonCodec.cs
Reed Solomon符号化実装クラス(別版)
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace iPentec.ReedSolomonCodec
{
public class iPentecReedSolomonCodec
{
private int NPAR = 0;
private int MAXDEG = 0;
private int FDimension = 0;
private iPentecGalois GT = null;
private int FNErasures = 0;
private int[] ErasureLocs;
private int[] ErrorLocs;
private int NErrors = 0;
private int[] Lambda;
// [MAXDEG];
private int[] Omega;
public int[] pBytes;
public int[] synBytes;
public int[] genPoly;
public int Dimension
{
get
{
return FDimension;
}
set
{
// [MAXDEG];
if (FDimension != value) {
FDimension = value;
}
}
}
//public TiPentecReedSolomonCodec(Component Aowner) : base(Aowner)
public iPentecReedSolomonCodec()
{
FDimension = 8;
GT = new iPentecGalois();
}
/*
protected override void Dispose()
{
}
*/
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// 初期化
public void initialize_ecc(int ParityLength)
{
NPAR = ParityLength;
MAXDEG = NPAR * 2;
// * Initialize the galois field arithmetic tables */
GT.InitTables(FDimension);
// init_galois_tables();
// * Compute the encoder generator polynomial */
genPoly = new int[MAXDEG * 2];
compute_genpoly(NPAR, out genPoly);
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// Decode
//デコード関数
//デコード関数を呼び出したのち check_syndrome でシンドロームを求める
public void decode_data(ref char[] data, int nbytes)
{
int i;
int j;
int sum;
synBytes = new int[MAXDEG];
for (j = 0; j < NPAR; j++) {
// for j:=0 to 8-1 do begin
sum = 0;
for (i = 0; i < nbytes; i++) {
sum = ((byte)data[i]) ^ GT.gmult(GT.gexp[j + 1], sum);
}
synBytes[j] = sum;
}
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// Check Syndrome
// シンドローム算出
// Decodeを事前に呼び出す必要がある
public int check_syndrome()
{
int result;
int i;
int nz;
nz = 0;
for (i = 0; i < NPAR; i++) {
if (synBytes[i] != 0) {
nz = 1;
}
}
result = nz;
return result;
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// Encode
// エンコード関数 コードワード版
public void encode_data(char[] msg, int nbytes, out char[] dst)
{
int i;
int dbyte;
int j;
int[] LFSR;
LFSR = new int[NPAR + 1];
for (i = 0; i <= NPAR; i++) {
LFSR[i] = 0;
}
for (i = 0; i < nbytes; i++) {
dbyte = ((byte)msg[i]) ^ LFSR[NPAR - 1];
for (j = NPAR - 1; j >= 1; j--) {
LFSR[j] = LFSR[j - 1] ^ GT.gmult(genPoly[j], dbyte);
}
LFSR[0] = GT.gmult(genPoly[0], dbyte);
}
pBytes = new int[NPAR];
for (i = 0; i < NPAR; i++) {
pBytes[i] = LFSR[i];
}
build_codeword(msg, nbytes, out dst);
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// Encode Parity
// エンコード関数 パリティ分離版
public void encode_data_parity(char[] msg, int nbytes, out char[] Parity)
{
int i;
int dbyte;
int j;
int[] LFSR;
LFSR = new int[NPAR + 1];
for (i = 0; i <= NPAR; i++) {
LFSR[i] = 0;
}
for (i = 0; i < nbytes; i++) {
dbyte = ((byte)msg[i]) ^ LFSR[NPAR - 1];
for (j = NPAR - 1; j >= 1; j--) {
LFSR[j] = LFSR[j - 1] ^ GT.gmult(genPoly[j], dbyte);
}
LFSR[0] = GT.gmult(genPoly[0], dbyte);
}
Parity = new char[NPAR];
pBytes = new int[NPAR];
for (i = 0; i < NPAR; i++) {
Parity[i] = ((char)LFSR[i]);
}
// for i := 0 to NPAR-1 do Pareity := pBytes[i];
}
public short crc_ccitt(string msg, int leng)
{
short result = 0;
return result;
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// correct errors erasures
// 誤り訂正処理
public int correct_errors_erasures(ref char[] codeword, int csize, int nerasures, int[] erasures)
{
int result;
int r;
int i;
int j;
byte err;
int num;
int denom;
FNErasures = nerasures;
for (i = 0; i < nerasures; i++) {
ErasureLocs[i] = erasures[i];
}
Modified_Berlekamp_Massey();
Find_Roots();
if ((NErrors <= NPAR) && (NErrors > 0)) {
// /* first check for illegal error locs */
for (r = 0; r < NErrors; r++) {
if (ErrorLocs[r] >= csize) {
result = 0;
return result;
}
}
for (r = 0; r < NErrors; r++) {
i = ErrorLocs[r];
// /* evaluate Omega at alpha^(-i) */
num = 0;
for (j = 0; j < MAXDEG; j++) {
num = num ^ GT.gmult(Omega[j], GT.gexp[((255 - i) * j) % 255]);
}
// /* evaluate Lambda' (derivative) at alpha^(-i) ; all odd powers disappear */
denom = 0;
j = 1;
while ((j < MAXDEG)) {
denom = denom ^ GT.gmult(Lambda[j], GT.gexp[((255 - i) * (j - 1)) % 255]);
j += 2;
}
err = (byte)GT.gmult(num, GT.ginv(denom));
//codeword[csize - i - 1] = ((char)((byte)codeword[csize - i - 1]) ^ err);
byte b = (byte)codeword[csize - i - 1];
b = (byte)(b ^ err);
codeword[csize - i - 1] = (char)b;
}
result = 1;
return result;
}
else {
result = 0;
}
return result;
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
// Compute GenPoly
// 誤り訂正関数作成
// 誤り訂正関数はエンコード時のみ使用される
protected void compute_genpoly(int nbytes, out int[] genpoly)
{
genpoly = new int[0];
int i;
int[] tp = new int[255 + 1];
int[] tp1 = new int[255 + 1];
// /* multiply (x + a^n) for n = 1 to nbytes */
zero_poly(out tp1);
tp1[0] = 1;
for (i = 1; i <= nbytes; i++) {
zero_poly(out tp);
tp[0] = GT.gexp[i];
// * set up x+a^n */
tp[1] = 1;
mult_polys(out genpoly, tp, tp1);
copy_poly(out tp1, genpoly);
}
// QRコード用 誤り訂正多項式
// for i := 0 to nbytes do begin
// genpoly[i]:=GT.gexp[GT.glog[genpoly[i]] - (nbytes-i)];
// end;
}
protected void build_codeword(char[] msg, int nbytes, out char[] dst)
{
int i;
dst = new char[nbytes];
for (i = 0; i < nbytes; i++) {
dst[i] = msg[i];
}
for (i = 0; i < NPAR; i++) {
dst[i + nbytes] = ((char)pBytes[NPAR - 1 - i]);
}
}
protected void Modified_Berlekamp_Massey()
{
int n;
int l;
int l2;
int k;
int d;
int i;
int[] psi;
int[] psi2;
int[] DA;
int[] gamma;
psi = new int[MAXDEG];
psi2 = new int[MAXDEG];
DA = new int[MAXDEG];
gamma = new int[MAXDEG];
// /* initialize Gamma, the erasure locator polynomial */
init_gamma(out gamma);
// /* initialize to z */
copy_poly(out DA, gamma);
mul_z_poly(ref DA);
copy_poly(out psi, gamma);
k = -1;
l = FNErasures;
for (n = FNErasures; n < NPAR; n++) {
// for n:=FNErasures to 8-1 do begin
d = compute_discrepancy(psi, synBytes, l, n);
if (d != 0) {
// /* psi2 = psi - d*D */
for (i = 0; i < MAXDEG; i++) {
psi2[i] = psi[i] ^ GT.gmult(d, DA[i]);
}
if (l < (n - k)) {
l2 = n - k;
k = n - l;
// /* D = scale_poly(ginv(d), psi); */
for (i = 0; i < MAXDEG; i++) {
DA[i] = GT.gmult(psi[i], GT.ginv(d));
}
l = l2;
}
// /* psi = psi2 */
for (i = 0; i < MAXDEG; i++) {
psi[i] = psi2[i];
}
}
mul_z_poly(ref DA);
}
// λの初期化
Lambda = new int[MAXDEG];
for (i = 0; i < MAXDEG; i++) {
Lambda[i] = psi[i];
}
compute_modified_omega();
}
// * gamma = product (1-z*a^Ij) for erasure locs Ij */
protected void init_gamma(out int[] gamma)
{
int e;
int[] tmp;
tmp = new int[MAXDEG];
zero_poly(out gamma);
zero_poly(out tmp);
gamma[0] = 1;
for (e = 0; e < FNErasures; e++) {
copy_poly(out tmp, gamma);
scale_poly(GT.gexp[ErasureLocs[e]], out tmp);
mul_z_poly(ref tmp);
add_polys(out gamma, tmp);
}
}
protected int compute_discrepancy(int[] lambda, int[] S, int L, int n)
{
int result;
int i;
int sum;
sum = 0;
for (i = 0; i <= L; i++) {
sum = sum ^ GT.gmult(lambda[i], S[n - i]);
}
result = sum;
return result;
}
protected void compute_modified_omega()
{
int i;
int[] product;
product = new int[MAXDEG * 2];
mult_polys(out product, Lambda, synBytes);
// ωの初期化 (0でクリア後 productで初期値を代入)
Omega = new int[MAXDEG];
zero_poly(out Omega);
for (i = 0; i < NPAR; i++) {
Omega[i] = product[i];
}
//this.Finalize(product);
}
protected void Find_Roots()
{
int sum;
int r;
int k;
NErrors = 0;
for (r = 1; r < 256; r++) {
sum = 0;
// /* evaluate lambda at r */
for (k = 0; k < NPAR + 1; k++) {
sum = sum ^ GT.gmult(GT.gexp[(k * r) % 255], Lambda[k]);
}
if ((sum == 0)) {
ErrorLocs[NErrors] = 255 - r;
NErrors++;
}
}
}
public void zero_poly(out int[] poly)
{
int i;
poly = new int[MAXDEG];
for (i = 0; i < MAXDEG; i++) {
poly[i] = 0;
}
}
public void copy_poly(out int[] dst, int[] src)
{
int i;
dst = new int[MAXDEG];
for (i = 0; i < MAXDEG; i++) {
dst[i] = src[i];
}
}
public void add_polys(out int[] dst, int[] src)
{
int i;
dst = new int[MAXDEG];
for (i = 0; i < MAXDEG; i++) {
dst[i] = dst[i] ^ src[i];
}
}
// /* multiply by z, i.e., shift right by 1 */
protected void mul_z_poly(ref int[] src)
{
int i;
src = new int[MAXDEG - 1];
for (i = MAXDEG - 1; i >= 1; i--) {
src[i] = src[i - 1];
}
src[0] = 0;
}
public void scale_poly(int k, out int[] poly)
{
int i;
poly = new int[MAXDEG];
for (i = 0; i < MAXDEG; i++) {
poly[i] = GT.gmult(k, poly[i]);
}
}
public void mult_polys(out int[] dst, int[] p1, int[] p2)
{
int i;
int j;
int[] tmp1;
tmp1 = new int[MAXDEG * 2];
dst = new int[MAXDEG * 2];
for (i = 0; i < (MAXDEG * 2); i++) {
dst[i] = 0;
}
for (i = 0; i < MAXDEG; i++) {
for (j = MAXDEG; j < (MAXDEG * 2); j++) {
tmp1[j] = 0;
}
// /* scale tmp1 by p1[i] */
for (j = 0; j < MAXDEG; j++) {
tmp1[j] = GT.gmult(p2[j], p1[i]);
}
// /* and mult (shift) tmp1 right by i */
j = (MAXDEG * 2) - 1;
while ((j >= i)) {
tmp1[j] = tmp1[j - i];
j -= 1;
}
for (j = 0; j < i; j++) {
tmp1[j] = 0;
}
// /* add into partial product */
for (j = 0; j < (MAXDEG * 2); j++) {
dst[j] = dst[j] ^ tmp1[j];
}
}
}
} // end TiPentecReedSolomonCodec
}
著者
iPentec Document 編集部
iPentec Document 編集部です。
快適な生活のための情報、価値のある体験やレビューを提供します。
最終更新日: 2024-01-24
改訂日: 2022-06-06
作成日: 2010-01-25
