| /* rs.c */ |
| /* This program is an encoder/decoder for Reed-Solomon codes. Encoding is in |
| systematic form, decoding via the Berlekamp iterative algorithm. |
| In the present form , the constants mm, nn, tt, and kk=nn-2tt must be |
| specified (the double letters are used simply to avoid clashes with |
| other n,k,t used in other programs into which this was incorporated!) |
| Also, the irreducible polynomial used to generate GF(2**mm) must also be |
| entered -- these can be found in Lin and Costello, and also Clark and Cain. |
| |
| The representation of the elements of GF(2**m) is either in index form, |
| where the number is the power of the primitive element alpha, which is |
| convenient for multiplication (add the powers modulo 2**m-1) or in |
| polynomial form, where the bits represent the coefficients of the |
| polynomial representation of the number, which is the most convenient form |
| for addition. The two forms are swapped between via lookup tables. |
| This leads to fairly messy looking expressions, but unfortunately, there |
| is no easy alternative when working with Galois arithmetic. |
| |
| The code is not written in the most elegant way, but to the best |
| of my knowledge, (no absolute guarantees!), it works. |
| However, when including it into a simulation program, you may want to do |
| some conversion of global variables (used here because I am lazy!) to |
| local variables where appropriate, and passing parameters (eg array |
| addresses) to the functions may be a sensible move to reduce the number |
| of global variables and thus decrease the chance of a bug being introduced. |
| |
| This program does not handle erasures at present, but should not be hard |
| to adapt to do this, as it is just an adjustment to the Berlekamp-Massey |
| algorithm. It also does not attempt to decode past the BCH bound -- see |
| Blahut "Theory and practice of error control codes" for how to do this. |
| |
| Simon Rockliff, University of Adelaide 21/9/89 |
| |
| 26/6/91 Slight modifications to remove a compiler dependent bug which hadn't |
| previously surfaced. A few extra comments added for clarity. |
| Appears to all work fine, ready for posting to net! |
| |
| Notice |
| -------- |
| This program may be freely modified and/or given to whoever wants it. |
| A condition of such distribution is that the author's contribution be |
| acknowledged by his name being left in the comments heading the program, |
| however no responsibility is accepted for any financial or other loss which |
| may result from some unforseen errors or malfunctioning of the program |
| during use. |
| Simon Rockliff, 26th June 1991 |
| */ |
| |
| #include <math.h> |
| #include <stdio.h> |
| #include <stdlib.h> |
| #define mm 8 /* RS code over GF(2**4) - change to suit */ |
| #define nn 255 /* nn=2**mm -1 length of codeword */ |
| #define tt 8 /* number of errors that can be corrected */ |
| #define kk 239 /* kk = nn-2*tt */ |
| |
| static int pp [mm+1] = { 1, 0, 1, 1, 1, 0, 0, 0, 1} ; /* specify irreducible polynomial coeffts */ |
| static int alpha_to [nn+1], index_of [nn+1], gg [nn-kk+1] ; |
| static int recd [nn], data [kk], bb [nn-kk] ; |
| static inited = 0; |
| |
| static void generate_gf() |
| /* generate GF(2**mm) from the irreducible polynomial p(X) in pp[0]..pp[mm] |
| lookup tables: index->polynomial form alpha_to[] contains j=alpha**i; |
| polynomial form -> index form index_of[j=alpha**i] = i |
| alpha=2 is the primitive element of GF(2**mm) |
| */ |
| { |
| register int i, mask ; |
| |
| mask = 1 ; |
| alpha_to[mm] = 0 ; |
| for (i=0; i<mm; i++) |
| { alpha_to[i] = mask ; |
| index_of[alpha_to[i]] = i ; |
| if (pp[i]!=0) |
| alpha_to[mm] ^= mask ; |
| mask <<= 1 ; |
| } |
| index_of[alpha_to[mm]] = mm ; |
| mask >>= 1 ; |
| for (i=mm+1; i<nn; i++) |
| { if (alpha_to[i-1] >= mask) |
| alpha_to[i] = alpha_to[mm] ^ ((alpha_to[i-1]^mask)<<1) ; |
| else alpha_to[i] = alpha_to[i-1]<<1 ; |
| index_of[alpha_to[i]] = i ; |
| } |
| index_of[0] = -1 ; |
| } |
| |
| |
| static void gen_poly() |
| /* Obtain the generator polynomial of the tt-error correcting, length |
| nn=(2**mm -1) Reed Solomon code from the product of (X+alpha**i), i=1..2*tt |
| */ |
| { |
| register int i,j ; |
| |
| gg[0] = 2 ; /* primitive element alpha = 2 for GF(2**mm) */ |
| gg[1] = 1 ; /* g(x) = (X+alpha) initially */ |
| for (i=2; i<=nn-kk; i++) |
| { gg[i] = 1 ; |
| for (j=i-1; j>0; j--) |
| if (gg[j] != 0) gg[j] = gg[j-1]^ alpha_to[(index_of[gg[j]]+i)%nn] ; |
| else gg[j] = gg[j-1] ; |
| gg[0] = alpha_to[(index_of[gg[0]]+i)%nn] ; /* gg[0] can never be zero */ |
| } |
| /* convert gg[] to index form for quicker encoding */ |
| for (i=0; i<=nn-kk; i++) gg[i] = index_of[gg[i]] ; |
| } |
| |
| |
| static void encode_rs() |
| /* take the string of symbols in data[i], i=0..(k-1) and encode systematically |
| to produce 2*tt parity symbols in bb[0]..bb[2*tt-1] |
| data[] is input and bb[] is output in polynomial form. |
| Encoding is done by using a feedback shift register with appropriate |
| connections specified by the elements of gg[], which was generated above. |
| Codeword is c(X) = data(X)*X**(nn-kk)+ b(X) */ |
| { |
| register int i,j ; |
| int feedback ; |
| |
| for (i=0; i<nn-kk; i++) bb[i] = 0 ; |
| for (i=kk-1; i>=0; i--) |
| { feedback = index_of[data[i]^bb[nn-kk-1]] ; |
| if (feedback != -1) |
| { for (j=nn-kk-1; j>0; j--) |
| if (gg[j] != -1) |
| bb[j] = bb[j-1]^alpha_to[(gg[j]+feedback)%nn] ; |
| else |
| bb[j] = bb[j-1] ; |
| bb[0] = alpha_to[(gg[0]+feedback)%nn] ; |
| } |
| else |
| { for (j=nn-kk-1; j>0; j--) |
| bb[j] = bb[j-1] ; |
| bb[0] = 0 ; |
| } ; |
| } ; |
| } ; |
| |
| |
| |
| static void decode_rs() |
| /* assume we have received bits grouped into mm-bit symbols in recd[i], |
| i=0..(nn-1), and recd[i] is index form (ie as powers of alpha). |
| We first compute the 2*tt syndromes by substituting alpha**i into rec(X) and |
| evaluating, storing the syndromes in s[i], i=1..2tt (leave s[0] zero) . |
| Then we use the Berlekamp iteration to find the error location polynomial |
| elp[i]. If the degree of the elp is >tt, we cannot correct all the errors |
| and hence just put out the information symbols uncorrected. If the degree of |
| elp is <=tt, we substitute alpha**i , i=1..n into the elp to get the roots, |
| hence the inverse roots, the error location numbers. If the number of errors |
| located does not equal the degree of the elp, we have more than tt errors |
| and cannot correct them. Otherwise, we then solve for the error value at |
| the error location and correct the error. The procedure is that found in |
| Lin and Costello. For the cases where the number of errors is known to be too |
| large to correct, the information symbols as received are output (the |
| advantage of systematic encoding is that hopefully some of the information |
| symbols will be okay and that if we are in luck, the errors are in the |
| parity part of the transmitted codeword). Of course, these insoluble cases |
| can be returned as error flags to the calling routine if desired. */ |
| { |
| register int i,j,u,q ; |
| int elp[nn-kk+2][nn-kk], d[nn-kk+2], l[nn-kk+2], u_lu[nn-kk+2], s[nn-kk+1] ; |
| int count=0, syn_error=0, root[tt], loc[tt], z[tt+1], err[nn], reg[tt+1] ; |
| |
| /* first form the syndromes */ |
| for (i=1; i<=nn-kk; i++) |
| { s[i] = 0 ; |
| for (j=0; j<nn; j++) |
| if (recd[j]!=-1) |
| s[i] ^= alpha_to[(recd[j]+i*j)%nn] ; /* recd[j] in index form */ |
| /* convert syndrome from polynomial form to index form */ |
| if (s[i]!=0) syn_error=1 ; /* set flag if non-zero syndrome => error */ |
| s[i] = index_of[s[i]] ; |
| } ; |
| |
| if (syn_error) /* if errors, try and correct */ |
| { |
| /* printf("RS: errors detected\n"); */ |
| /* 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. |
| */ |
| /* initialise 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<nn-kk; 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 ; |
| |
| do |
| { |
| 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] = index_of[elp[u][i]] ; |
| } |
| } |
| else |
| /* search for words with greatest u_lu[q] for which d[q]!=0 */ |
| { q = u-1 ; |
| while ((d[q]==-1) && (q>0)) q-- ; |
| /* have found first non-zero d[q] */ |
| if (q>0) |
| { j=q ; |
| do |
| { j-- ; |
| if ((d[j]!=-1) && (u_lu[q]<u_lu[j])) |
| q = j ; |
| }while (j>0) ; |
| } ; |
| |
| /* 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] = l[q]+u-q ; |
| |
| /* form new elp(x) */ |
| for (i=0; i<nn-kk; i++) elp[u+1][i] = 0 ; |
| for (i=0; i<=l[q]; i++) |
| if (elp[q][i]!=-1) |
| elp[u+1][i+u-q] = alpha_to[(d[u]+nn-d[q]+elp[q][i])%nn] ; |
| for (i=0; i<=l[u]; i++) |
| { elp[u+1][i] ^= elp[u][i] ; |
| elp[u][i] = index_of[elp[u][i]] ; /*convert old elp value to index*/ |
| } |
| } |
| u_lu[u+1] = u-l[u+1] ; |
| |
| /* form (u+1)th discrepancy */ |
| if (u<nn-kk) /* no discrepancy computed on last iteration */ |
| { |
| if (s[u+1]!=-1) |
| d[u+1] = alpha_to[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] ^= alpha_to[(s[u+1-i]+index_of[elp[u+1][i]])%nn] ; |
| d[u+1] = index_of[d[u+1]] ; /* put d[u+1] into index form */ |
| } |
| } while ((u<nn-kk) && (l[u+1]<=tt)) ; |
| |
| u++ ; |
| if (l[u]<=tt) /* can correct error */ |
| { |
| /* put elp into index form */ |
| for (i=0; i<=l[u]; i++) elp[u][i] = index_of[elp[u][i]] ; |
| |
| /* find roots of the error location polynomial */ |
| for (i=1; i<=l[u]; i++) |
| reg[i] = elp[u][i] ; |
| count = 0 ; |
| for (i=1; i<=nn; i++) |
| { q = 1 ; |
| for (j=1; j<=l[u]; j++) |
| if (reg[j]!=-1) |
| { reg[j] = (reg[j]+j)%nn ; |
| q ^= alpha_to[reg[j]] ; |
| } ; |
| if (!q) /* store root and error location number indices */ |
| { root[count] = i; |
| loc[count] = nn-i ; |
| count++ ; |
| }; |
| } ; |
| if (count==l[u]) /* no. roots = degree of elp hence <= tt errors */ |
| { |
| /* 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] = alpha_to[s[i]] ^ alpha_to[elp[u][i]] ; |
| else if ((s[i]!=-1) && (elp[u][i]==-1)) |
| z[i] = alpha_to[s[i]] ; |
| else if ((s[i]==-1) && (elp[u][i]!=-1)) |
| z[i] = alpha_to[elp[u][i]] ; |
| else |
| z[i] = 0 ; |
| for (j=1; j<i; j++) |
| if ((s[j]!=-1) && (elp[u][i-j]!=-1)) |
| z[i] ^= alpha_to[(elp[u][i-j] + s[j])%nn] ; |
| z[i] = index_of[z[i]] ; /* put into index form */ |
| } ; |
| |
| /* evaluate errors at locations given by error location numbers loc[i] */ |
| for (i=0; i<nn; i++) |
| { err[i] = 0 ; |
| if (recd[i]!=-1) /* convert recd[] to polynomial form */ |
| recd[i] = alpha_to[recd[i]] ; |
| else recd[i] = 0 ; |
| } |
| for (i=0; i<l[u]; 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]] ^= alpha_to[(z[j]+j*root[i])%nn] ; |
| if (err[loc[i]]!=0) |
| { err[loc[i]] = index_of[err[loc[i]]] ; |
| q = 0 ; /* form denominator of error term */ |
| for (j=0; j<l[u]; j++) |
| if (j!=i) |
| q += index_of[1^alpha_to[(loc[j]+root[i])%nn]] ; |
| q = q % nn ; |
| err[loc[i]] = alpha_to[(err[loc[i]]-q+nn)%nn] ; |
| recd[loc[i]] ^= err[loc[i]] ; /*recd[i] must be in polynomial form */ |
| } |
| } |
| } |
| else /* no. roots != degree of elp => >tt errors and cannot solve */ |
| for (i=0; i<nn; i++) /* could return error flag if desired */ |
| if (recd[i]!=-1) /* convert recd[] to polynomial form */ |
| recd[i] = alpha_to[recd[i]] ; |
| else recd[i] = 0 ; /* just output received codeword as is */ |
| } |
| else /* elp has degree has degree >tt hence cannot solve */ |
| for (i=0; i<nn; i++) /* could return error flag if desired */ |
| if (recd[i]!=-1) /* convert recd[] to polynomial form */ |
| recd[i] = alpha_to[recd[i]] ; |
| else recd[i] = 0 ; /* just output received codeword as is */ |
| } |
| else /* no non-zero syndromes => no errors: output received codeword */ |
| for (i=0; i<nn; i++) |
| if (recd[i]!=-1) /* convert recd[] to polynomial form */ |
| recd[i] = alpha_to[recd[i]] ; |
| else recd[i] = 0 ; |
| } |
| |
| |
| void rsdec_204(unsigned char* data_out, unsigned char* data_in) |
| { |
| int i; |
| |
| if (!inited) { |
| /* generate the Galois Field GF(2**mm) */ |
| generate_gf(); |
| /* compute the generator polynomial for this RS code */ |
| gen_poly(); |
| |
| inited = 1; |
| } |
| |
| /* put the transmitted codeword, made up of data plus parity, in recd[] */ |
| |
| /* parity */ |
| for (i=0; i<204-188; ++i) { |
| recd[i] = data_in[188 + i]; |
| } |
| /* zeroes */ |
| for (i=0; i<255-204; ++i) { |
| recd[204-188 + i] = 0; |
| } |
| /* data */ |
| for (i=0; i<188; ++i) { |
| recd[255-188 + i] = data_in[i]; |
| } |
| |
| for (i=0; i<nn; i++) |
| recd[i] = index_of[recd[i]] ; /* put recd[i] into index form */ |
| |
| /* decode recv[] */ |
| decode_rs() ; /* recd[] is returned in polynomial form */ |
| |
| for (i=0; i<188; ++i) { |
| data_out [i] = recd[255-188 + i]; |
| } |
| } |
| |
| void rsenc_204(unsigned char* data_out, unsigned char* data_in) |
| { |
| int i; |
| |
| if (!inited) { |
| /* generate the Galois Field GF(2**mm) */ |
| generate_gf(); |
| /* compute the generator polynomial for this RS code */ |
| gen_poly(); |
| |
| inited = 1; |
| } |
| |
| /* zeroes */ |
| for (i=0; i<255-204; ++i) { |
| data[i] = 0; |
| } |
| /* data */ |
| for (i=0; i<188; ++i) { |
| data[255-204 + i] = data_in[i]; |
| } |
| |
| encode_rs(); |
| |
| for (i=0; i<188; ++i) { |
| data_out[i] = data_in[i]; |
| } |
| for (i=0; i<204-188; ++i) { |
| data_out[188+i] = bb[i]; |
| } |
| |
| } |
| |
| int main(void) { |
| unsigned char rs_in[204], rs_out[204]; |
| int i, j, k; |
| |
| #ifdef SMALL_PROBLEM_SIZE |
| #define LENGTH 15000 |
| #else |
| #define LENGTH 150000 |
| #endif |
| |
| #ifdef __MINGW32__ |
| #define random() rand() |
| #endif |
| |
| for (i=0; i<LENGTH; ++i) { |
| /* Generate random data */ |
| for (j=0; j<188; ++j) { |
| rs_in[j] = (random() & 0xFF); |
| } |
| rsenc_204(rs_out, rs_in); |
| /* Number of errors to insert */ |
| k = random() & 0x7F; |
| |
| for (j=0; j<k; ++j) { |
| rs_out[random() % 204] = (random() & 0xFF); |
| } |
| |
| rsdec_204(rs_in, rs_out); |
| } |
| return 0; |
| } |
