/*******************************************************************************
*                                                                              *
*  File:       createTable.c                            Revision:  1.0         *
*                                                                              *
*  Contents:   a simple program that creates RSE encoding/decoding tables that *
*              can then be used for a Java-based RSE encoder/decoder           *
*                                                                              *
*  Creation:   05.02.1998                      Last Modification:  13.02.1998  *
*                                                                              *
*  Platform:   SGI R4000 Indigo running IRIX 5.3                               *
*                                                                              *
*  Program:    ANSI C                                                          *
*                                                                              *
*  Author:     Andreas Rozek                                                   *
*              Stuttgart University Computer Center                            *
*              Communication Systems & BelWü Development                       *
*              Allmandring 3a                                                  *
*            D-70550 Stuttgart                                                 *
*              Germany                                                         *
*                                                                              *
*    Phone:    ++49 (711) 685-4514                                             *
*    Fax:      ++49 (711) 678-8363                                             *
*    EMail:    Andreas.Rozek@RUS.Uni-Stuttgart.De                              *
*                                                                              *
*  Comments:   the program should be invoked as follows:                       *
*                                                                              *
*                createTable <filename> <symbolsize> <blocksize> <paritycount> *
*                                                                              *
*              with the following meanings for the command line arguments:     *
*                                                                              *
*                <filename>     name of the file to write the table into       *
*                <symbolsize>   size of code symbols in bits                   *
*                <blocksize>    number of symbols per transmission block       *
*                               (must be <= (2**symbolsize)-1)                 *
*                               Note:  this program is able  to produce tables *
*                               for  "shortened" RSE codes  (i.e.,  codes with *
*                               block sizes < (2**symbolsize)-1)               *
*                <paritycount>  number of parity symbols per block  (that many *
*                               symbol losses  may be repared by the resulting *
*                               RSE code)                                      *
*                                                                              *
*              The final RSE encoding/decoding is based on code written by     *
*                                                                              *
*                Phil Karn (karn@ka9q.ampr.org)                                *
*                                                                              *
*              and was only changed to become  a little bit more flexible than *
*              the original                                                    *
*                                                                              *
*******************************************************************************/

#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>

/*******************************************************************************
*                                                                              *
*                             Function Prototypes                              *
*                                                                              *
*******************************************************************************/

  void showUsage (char *PrgName);

  int  Cases (int StartIndex, int SymbolCount, int Size);

  void writeByte  (FILE *OutFile, unsigned int Value);
  void writeShort (FILE *OutFile, unsigned int Value);
  void writeInt   (FILE *OutFile, unsigned int Value);

  void initCodec (int SymbolSize);
  void definePrimitivePolynomial (int Count, ...);

  void generateGaloisField (void);
  void generatePolynomial  (void);

  int  encode (int DataSymbol[], int ParitySymbol[]);
  int  decode (int DataSymbol[], int LossList[], int ParityCount);

/*******************************************************************************
*                                                                              *
*                               Global Variables                               *
*                                                                              *
*******************************************************************************/

  int   SymbolSize;                           /* size of code symbols in bits */
  int   BlockSize;                  /* size of transmission blocks in symbols */
  int   DataCount;           /* number of data symbols per transmission block */
  int   ParityCount;       /* number of parity symbols per transmission block */

  int   FullBlockSize;        /* full block size (2**SymbolSize-1) in symbols */
  int   FullDataCount;           /* # data symbols for "unshortened" RSE code */

  int   PrimitivePolynomial[17];                      /* primitive polynomial */

/*******************************************************************************
*                                                                              *
* main                                                            main program *
*                                                                              *
*******************************************************************************/

  int main (int argc, char *argv[]) {
    char *FileName;              /* name of the file to write the tables into */

    int   ParityBits;            /* size of encoding table elements (in bits) */
    int   DataBits;             /* size of encoding table "address" (in bits) */

    int   CaseCount;   /* # of failure cases that can be handled by this code */
    int   DataCases, ParityCases;                           /* loop variables */

    int   EncoderLength, DecoderLength; /* # of encoder/decoder table entries */
    int  *EncoderTable, *DecoderTable;       /* actual encoder/decoder tables */

    int   SymbolMask;              /* to mask out contents of a single symbol */
    int  *DataSymbol, *ParitySymbol;     /* in/out data of internal RSE codec */
    int   Parity;                           /* a single encoder table content */

    int  *LossList, *ParityList;             /* list of loss/parity positions */
    int   LossCount;                                /* number of symbols lost */
    int   LossBase, LossFactor;/* numbers used to calculate a loss case index */
    int   LossCase, UseCase;                                          /* dto. */
    int   Position;              /* symbol position when traversing the block */

    int   Result;          /* result of RSE ecoding with given parity symbols */

    FILE *OutFile;                           /* file to write the tables into */

    int   i,j,k,l,m,n;                                      /* loop variables */


    setbuf(stdout, NULL);                      /* for debugging purposes only */

    printf("%s - RSE code table generator\n", argv[0]);

  /**** fetch command line arguments ****/

    if (argc != 5) {
      showUsage(argv[0]);
      exit(1);
    };

    FileName    = argv[1];
    SymbolSize  = atoi(argv[2]);
    BlockSize   = atoi(argv[3]);
    ParityCount = atoi(argv[4]);

  /**** now check the given numbers (regardless of any later restrictions) ****/

    if ((SymbolSize < 2) || (SymbolSize > 16)) {
      printf("  illegal SymbolSize \"%s\"\n", argv[2]);
      printf("  please specify a numerical value in the range 2...16\n");
      exit(1);
    };

    if ((BlockSize < 2) || (BlockSize > 65535)) {
      printf("  illegal BlockSize \"%s\"\n", argv[3]);
      printf("  please specify a numerical value in the range 2...65535\n");
      exit(2);
    };

    if (BlockSize > (1 << SymbolSize)-1) {
      printf("  BlockSize %u is too large for symbols of size %u\n", BlockSize, SymbolSize);
      printf("  please specify a numerical value in the range 2...2**SymbolSize-1\n");
      exit(3);
    };

    if ((ParityCount < 1) || (ParityCount > 65535-1)) {
      printf("  illegal ParityCount \"%s\"\n", argv[4]);
      printf("  please specify a numerical value in the range 1...65535-1\n");
      exit(4);
    };

    if (ParityCount >= BlockSize) {
      printf("  ParityCount %u is too large for blocks of size %u\n", ParityCount, BlockSize);
      printf("  please specify a value in the range 1...BlockSize-1\n");
      exit(5);
    };

  /**** currently, there are several implementation limits ****/

    ParityBits = ParityCount*SymbolSize;
    if (ParityBits > 32) {
      printf("  ParityCount (%u) and SymbolSize (%u) are too large\n", ParityCount,SymbolSize);
      printf("  please ensure that ParityCount*SymbolSize <= 32\n");
      exit(10);
    };

    DataCount = BlockSize-ParityCount;
    DataBits  = DataCount*SymbolSize;
    if (DataBits > 20) {
      printf("  BlockSize-ParityCount (%u) and SymbolSize (%u) are too large\n", DataCount,SymbolSize);
      printf("  please ensure that (BlockSize-ParityCount)*SymbolSize <= 20\n");
      exit(11);
    };

    EncoderLength = 1 << DataBits;               /* size of encoder addresses */

  /**** determine the number of relevant failure cases ****/

    CaseCount = 0;
    for (i = 1; i <= (DataCount < ParityCount ? DataCount : ParityCount); i++) {
      DataCases = DataCount;
      for (j = 1; j < i; j++) DataCases   =   (DataCases * (DataCount-j))   / (j+1);

      ParityCases = ParityCount;
      for (j = 1; j < i; j++) ParityCases = (ParityCases * (ParityCount-j)) / (j+1);

      CaseCount += DataCases*ParityCases;
    };

    DecoderLength = (1 << DataBits) * CaseCount; /* size of decoder addresses */
    if (DecoderLength > (1 << 20)) {
      printf("  there are too many possible loss cases (%u)\n", CaseCount);
      printf("  try reducing BlockSize or ParityCount\n");
      exit(12);
    };

  /**** inform about important table parameters ****/

    printf("\n");
    printf("  creating RSE encoding/decoding tables with the following parameters:\n");
    printf("\n");
    printf("    Output File: \"%s\"\n", FileName);
    printf("\n");
    printf("    Symbol Size:          %6u Bits\n",      SymbolSize);
    printf("    Block Size:           %6u Symbols\n",   BlockSize);
    printf("    Data Symbols:         %6u per block\n", DataCount);
    printf("    Parity Symbols:       %6u per block\n", ParityCount);
    printf("\n");
    printf("  Encoding Table:\n");
    printf("\n");
    printf("    Number of Entries:    %6u\n",           EncoderLength);
    printf("    Size of each Entry:   %6u Bits\n",      ParityBits);
    printf("\n");
    printf("  Decoding Table:\n");
    printf("\n");
    printf("    Number of Loss Cases: %6u\n",           CaseCount);
    printf("    Number of Entries:    %6u\n",           DecoderLength);
    printf("    Size of each Entry:   %6u Bits\n",      ParityBits);

  /**** initialize RSE codec ****/

    FullBlockSize = (1 << SymbolSize)-1;                   /* 2**SymbolSize-1 */
    FullDataCount = FullBlockSize-ParityCount;

    initCodec(SymbolSize);

  /**** create RSE encoding table ****/

    printf("\n");
    printf("  creating RSE encoding table...\n");

    DataSymbol   = malloc(sizeof(int)*FullBlockSize); /* for both en/decoding */
    ParitySymbol = malloc(sizeof(int)*ParityCount);

    SymbolMask = (1 << SymbolSize)-1;    /* to mask symbols out of an integer */

    for (j = DataCount; j < FullBlockSize; j++) {
      DataSymbol[j] = 0;                  /* prepare for "shortened" RSE code */
    };

    EncoderTable = malloc(sizeof(int)*EncoderLength);
    for (i = 0; i < EncoderLength; i++) {

    /**** create data to be encoded (from EncoderTable index) ****/

      for (j = 0; j < DataCount; j++) {  /* scatter "i" into separate symbols */
        DataSymbol[DataCount-j-1] = (i >> j*SymbolSize) & SymbolMask;
      };

    /**** calculate RSE encoding ****/

      encode (DataSymbol, ParitySymbol);

    /**** produce encoder table entry ****/

      Parity = 0;
      for (j = 1; j <= ParityCount; j++) {/* gather parities into table entry */
        Parity = (Parity << SymbolSize) | (ParitySymbol[ParityCount-j] & SymbolMask);
      };

      EncoderTable[i] = Parity;
    };

    printf("  done\n");

  /**** create RSE decoding table ****/

    printf("\n");
    printf("  creating RSE decoding table...\n");

    LossList   = malloc(sizeof(int)*ParityCount);
    ParityList = malloc(sizeof(int)*ParityCount);

    DecoderTable = malloc(sizeof(int)*DecoderLength);
    for (i = 0; i < EncoderLength; i++) {
      for (j = 0; j < CaseCount; j++) {  /* handle every loss case separately */

      /**** calculate number of lost symbols ****/

        DataCases   = DataCount;                /* # cases with 1 lost symbol */
        ParityCases = ParityCount;           /* # choices for 1 parity symbol */
        LossFactor  = DataCases*ParityCases;
        LossBase    = 0;
        for (LossCount = 1; LossCount <= ParityCount; LossCount++) {
          if (j < LossBase+LossFactor) {    /* do we have "LossCount" losses? */
            break;
          } else {
            DataCases   =   (DataCases*(DataCount-LossCount))  /(LossCount+1);
            ParityCases = (ParityCases*(ParityCount-LossCount))/(LossCount+1);

            LossBase   += LossFactor;
            LossFactor *= DataCases*ParityCases;
          };
        };

        LossCase = j-LossBase;                        /* "normalize" LossCase */

      /**** calculate positions of lost symbols ****/

        UseCase  = LossCase % ParityCases;      /* isolate parity information */
        LossCase = LossCase / ParityCases;        /* isolate data information */

        Position = 0;                    /* start at the beginning of a block */
        for (k = 0; k < LossCount-1; k++) {/* all lost sym.s but the last one */
          LossBase = 0; /* step through "LossCase" and test every possibility */
          for (l = Position; l < DataCount-(LossCount-1)+k; l++) {
            LossFactor = Cases(l+1, LossCount-k-1, DataCount);
            if (LossCase < LossBase+LossFactor) {
              LossList[k] = l;       /* remember position of lost data symbol */
              Position = l+1;                      /* prepare for next symbol */
              LossCase -= LossBase;
              break;
            } else {
              LossBase += LossFactor;
            };
          };
        };
        LossList[LossCount-1] = Position+LossCase;/* the last case is trivial */

      /**** calculate positions of parity symbols to be used ****/

        LossCase = UseCase;     /* poor, but allows to "reuse" LossCase below */

        Position = 0;         /* start at the beginning of parity information */
        for (k = 0; k < LossCount-1; k++) {  /* every symbol but the last one */
          LossBase = 0; /* step through "LossCase" and test every possibility */
          for (l = Position; l < ParityCount-(LossCount-1)+k; l++) {
            LossFactor = Cases(l+1, LossCount-k-1, ParityCount);
            if (LossCase < LossBase+LossFactor) {
              ParityList[k] = l;   /* remember position of used parity symbol */
              Position = l+1;                      /* prepare for next symbol */
              LossCase -= LossBase;
              break;
            } else {
              LossBase += LossFactor;
            };
          };
        };
        ParityList[LossCount-1] = Position+LossCase;   /* the last is trivial */

      /**** create data to be decoded ****/

        m = SymbolSize*(DataCount-1);     /* for both data and parity symbols */
        for (k = 0, l = 0; k < DataCount; k++) {
          if (LossList[l] == k) {             /* has actual symbol been lost? */
            DataSymbol[k] = 0;                     /* prepare a "lost" symbol */
            if (l < LossCount-1) l++; /* if need be: look at next lost symbol */
          } else {
            DataSymbol[k] = (i >> m) & SymbolMask;       /* set normal symbol */
            m -= SymbolSize;            /* proceed to next data/parity symbol */
          };
        };

      /**** append parity information (and mark unused parity symbols) ****/

        for (k = 0, l = 0, n = LossCount; k < ParityCount; k++) {
          if (ParityList[l] == k) {    /* is actual parity symbol to be used? */
            DataSymbol[FullDataCount+k] = (i >> m) & SymbolMask;
            m -= SymbolSize;
            if (l < LossCount-1) l++; /* if need be: look at next parity sym. */
          } else {
            DataSymbol[FullDataCount+k] = 0; /* clear unused/lost parity sym. */
            LossList[n] = FullDataCount+k;         /* mark lost parity symbol */
            n++;
          };
        };

      /**** calculate RSE decoding ****/

        Result = decode (DataSymbol, LossList, ParityCount);

      /**** determine decoder table address ****/

        k = (j << DataBits) | i;

      /**** produce decoder table entry ****/

        DecoderTable[k] = 0;
        if (Result > -1) {                       /* successful reconstruction */
          for (l = 0; l < LossCount; l++) {
            DecoderTable[k] |= ((DataSymbol[LossList[l]] & SymbolMask) << (l * SymbolSize));
          };
        } else {          /* reconstruction failed - this should never happen */
          printf("  no reconstruction possible for i = %6u, j = %5u \n", i, j);
          printf("    data symbols:        ");
            for (l = 0; l < DataCount-1; l++) {
              printf("%3u, ", DataSymbol[l]);
            };
            printf("%3u\n", DataSymbol[DataCount-1]);
          printf("    parity symbols:      ");
            for (l = 0; l < LossCount-1; l++) {
              printf("%3u, ", DataSymbol[FullDataCount+l]);
            };
            printf("%3u\n", DataSymbol[FullDataCount+LossCount-1]);
          printf("    lost data symbols:   ");
            for (l = 0; l < LossCount-1; l++) {
              printf("%3u, ", LossList[l]);
            };
            printf("%3u\n", LossList[LossCount-1]);
          printf("    used parity symbols: ");
            for (l = 0; l < LossCount-1; l++) {
              printf("%3u, ", ParityList[l]);
            };
            printf("%3u\n", ParityList[LossCount-1]);

          exit (20);
        };
      };
    };

    printf("  done\n");

  /**** write tables onto file ****/

    printf("\n");
    printf("  writing tables onto file...\n");

    OutFile = fopen(FileName, "wb");
    if (OutFile == NULL) {
      printf("  unable to open file \"%s\" for writing\n", FileName);
      printf("  tables will not be written\n");
      exit(12);
    };

      writeByte(OutFile, SymbolSize);
      writeByte(OutFile, BlockSize);
      writeByte(OutFile, DataCount);
      writeByte(OutFile, ParityCount);

      printf("    encoder table: %6u ", EncoderLength);

      if        (ParityBits <=  8) {
        printf("bytes\n");
        for (i = 0; i < EncoderLength; i++) {
          writeByte (OutFile, EncoderTable[i]);
        };
      } else if (ParityBits <= 16) {
        printf("words\n");
        for (i = 0; i < EncoderLength; i++) {
          writeShort(OutFile, EncoderTable[i]);
        };
      } else {
        printf("quadwords\n");
        for (i = 0; i < EncoderLength; i++) {
          writeInt  (OutFile, EncoderTable[i]);
        };
      };

      printf("    decoder table: %6u ", DecoderLength);

      if        (ParityBits <=  8) {
        printf("bytes\n");
        for (i = 0; i < DecoderLength; i++) {
          writeByte (OutFile, DecoderTable[i]);
        };
      } else if (ParityBits <= 16) {
        printf("words\n");
        for (i = 0; i < DecoderLength; i++) {
          writeShort(OutFile, DecoderTable[i]);
        };
      } else {
        printf("quadwords\n");
        for (i = 0; i < DecoderLength; i++) {
          writeInt  (OutFile, DecoderTable[i]);
        };
      };

    fclose (OutFile);

    printf("  done\n");

  /**** clean up and leave... ****/

    printf("\n");
    printf("RSE table creation finished\n");

    exit(0);
  }

/*******************************************************************************
*                                                                              *
* showUsage                          prints some usage information on "stdout" *
*                                                                              *
*******************************************************************************/

 void showUsage (char *PrgName) {
    printf("  usage: %s <filename> <symbolsize> <blocksize> <paritycount>\n", PrgName);
    printf("\n");
    printf("  <filename>     name of the file to write the table into\n");
    printf("  <symbolsize>   size of code symbols in bits\n");
    printf("  <blocksize>    number of symbols per transmission block\n");
    printf("                 (must be <= (2**symbolsize)-1)\n");
    printf("                 Note:  this program is able  to produce tables\n");
    printf("                 for  \"shortened\" RSE codes  (i.e.,  codes with\n");
    printf("                 block sizes < (2**symbolsize)-1)\n");
    printf("  <paritycount>  number of parity symbols per block  (that many\n");
    printf("                 symbol losses may be repared  by the resulting\n");
    printf("                 RSE code)\n");
    printf("\n");
 }

/*******************************************************************************
*                                                                              *
* Cases                      calculates the number of remaining loss/use cases *
*                                                                              *
*******************************************************************************/

  int Cases (int StartIndex, int SymbolCount, int Size) {
    int CaseCount;                         /* is used to sum up several terms */
    int Index;                                               /* loop variable */

    if (SymbolCount == 1) {
      return Size-StartIndex;                   /* "StartIndex" begins with 0 */
    } else {
      SymbolCount--;
      CaseCount = 0;
        for (Index = StartIndex; Index < Size-SymbolCount; Index++) {
          CaseCount += Cases(Index+1, SymbolCount, Size);
        };
      return CaseCount;
    };
  }

/*******************************************************************************
*                                                                              *
* writeByte                                        write a byte into "OutFile" *
*                                                                              *
*******************************************************************************/

  void writeByte (FILE *OutFile, unsigned int Value) {
    putc(Value & 0xff, OutFile);
  }

/*******************************************************************************
*                                                                              *
* writeShort                                   writes two bytes into "OutFile" *
*                                                                              *
*******************************************************************************/

  void writeShort (FILE *OutFile, unsigned int Value) {
    putc(Value & 0xff, OutFile);
      Value >>= 8;
    putc(Value & 0xff, OutFile);
  }

/*******************************************************************************
*                                                                              *
* writeInt                                    writes four bytes into "OutFile" *
*                                                                              *
*******************************************************************************/

  void writeInt (FILE *OutFile, unsigned int Value) {
    putc(Value & 0xff, OutFile);
      Value >>= 8;
    putc(Value & 0xff, OutFile);
      Value >>= 8;
    putc(Value & 0xff, OutFile);
      Value >>= 8;
    putc(Value & 0xff, OutFile);
  }

/*******************************************************************************
*                                                                              *
*                             RSE encoder/decoder                              *
*                                                                              *
*******************************************************************************/

/*******************************************************************************
*                                                                              *
* definePrimitivePolynomial                       defines primitive polynomial *
*                                                                              *
*******************************************************************************/

  void definePrimitivePolynomial (int Count, ...) {
    va_list ParamList;
    int     i;

    va_start(ParamList, Count);
      for (i = 0; i < Count; i++) {
        PrimitivePolynomial[i] = va_arg(ParamList, int);
      };
    va_end(ParamList);
  }

/**** mappings of my variable names to those of Phil Karn ****/

#define MM SymbolSize
#define KK FullDataCount
#define NN FullBlockSize

/*******************************************************************************
*                                                                              *
*             RSE encoder/decoder by Phil Karn, slightly modified              *
*                                                                              *
*******************************************************************************/

#define B0 1 /* alpha exponent for the first root of the generator polynomial */

  int Alpha_to[1 << 16];           /* index->polynomial form conversion table */
  int Index_of[1 << 16];           /* Polynomial->index form conversion table */

/* No legal value in index form represents zero, so we need a special value for this purpose */

#define A0 (NN)

/* generator polynomial g(x); degree of g(x) = 2*TT; has roots @**B0, @**(B0+1), ... ,@^(B0+2*TT-1) */

  int GeneratorPolynomial[1 << 16];

/*******************************************************************************
*                                                                              *
* modnn           compute x % NN, where NN is 2**MM - 1, without a slow divide *
*                                                                              *
*******************************************************************************/

  static int modnn (int x) {
    while (x >= NN) {
      x -= NN;
      x = (x >> MM) + (x & NN);
    }
    return x;
  }

#define min(a,b) ((a) < (b) ? (a) : (b))

#define CLEAR(a,n) {        \
  int ci;                   \
  for(ci=(n)-1;ci >=0;ci--) \
    (a)[ci] = 0;            \
  }

#define COPY(a,b,n) {       \
  int ci;                   \
  for(ci=(n)-1;ci >=0;ci--) \
    (a)[ci] = (b)[ci];      \
  }
#define COPYDOWN(a,b,n) {   \
  int ci;                   \
  for(ci=(n)-1;ci >=0;ci--) \
    (a)[ci] = (b)[ci];      \
  }

/*******************************************************************************
*                                                                              *
* initCodec                                     initialize RSE encoder/decoder *
*                                                                              *
*******************************************************************************/

  void initCodec (int SymbolSize) {
    switch (SymbolSize) {
      case  2: definePrimitivePolynomial( 3, 1, 1, 1); break;
      case  3: definePrimitivePolynomial( 4, 1, 1, 0, 1); break;
      case  4: definePrimitivePolynomial( 5, 1, 1, 0, 0, 1); break;
      case  5: definePrimitivePolynomial( 6, 1, 0, 1, 0, 0, 1); break;
      case  6: definePrimitivePolynomial( 7, 1, 1, 0, 0, 0, 0, 1); break;
      case  7: definePrimitivePolynomial( 8, 1, 0, 0, 1, 0, 0, 0, 1); break;
      case  8: definePrimitivePolynomial( 9, 1, 0, 1, 1, 1, 0, 0, 0, 1); break;
      case  9: definePrimitivePolynomial(10, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1); break;
      case 10: definePrimitivePolynomial(11, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1); break;
      case 11: definePrimitivePolynomial(12, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1); break;
      case 12: definePrimitivePolynomial(13, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1); break;
      case 13: definePrimitivePolynomial(14, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1); break;
      case 14: definePrimitivePolynomial(15, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1); break;
      case 15: definePrimitivePolynomial(16, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1); break;
      case 16: definePrimitivePolynomial(17, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1); break;
    };

    generateGaloisField();
    generatePolynomial();
  }

/*******************************************************************************
*                                                                              *
* generate GF(2**m) from the irreducible polynomial p(X) in p[0]..p[m]         *
*  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**m)                                *
*  HARI's COMMENT: (4/13/94) alpha_to[] can be used as follows:                *
*       Let @ represent the primitive element commonly called "alpha" that     *
*  is the root of the primitive polynomial p(x). Then in GF(2^m), for any      *
*  0 <= i <= 2^m-2,                                                            *
*       @^i = a(0) + a(1) @ + a(2) @^2 + ... + a(m-1) @^(m-1)                  *
*  where the binary vector (a(0),a(1),a(2),...,a(m-1)) is the representation   *
*  of the integer "alpha_to[i]" with a(0) being the LSB and a(m-1) the MSB.    *
*  Thus for example the polynomial representation of @^5 would be given by the *
*  binary representation of the integer "alpha_to[5]".                         *
*                  Similarily, index_of[] can be used as follows:              *
*       As above, let @ represent the primitive element of GF(2^m) that is     *
*  the root of the primitive polynomial p(x). In order to find the power       *
*  of @ (alpha) that has the polynomial representation                         *
*       a(0) + a(1) @ + a(2) @^2 + ... + a(m-1) @^(m-1)                        *
*  we consider the integer "i" whose binary representation with a(0) being LSB *
*  and a(m-1) MSB is (a(0),a(1),...,a(m-1)) and locate the entry               *
*  "index_of[i]". Now, @^index_of[i] is that element whose polynomial          *
*   representation is (a(0),a(1),a(2),...,a(m-1)).                             *
*  NOTE:                                                                       *
*       The element alpha_to[2^m-1] = 0 always signifying that the             *
*  representation of "@^infinity" = 0 is (0,0,0,...,0).                        *
*       Similarily, the element index_of[0] = A0 always signifying             *
*  that the power of alpha which has the polynomial representation             *
*  (0,0,...,0) is "infinity".                                                  *
*                                                                              *
*******************************************************************************/

  void generateGaloisField (void) {
    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 PrimitivePolynomial[i] == 1 then term @^i occurs in poly-repr of @^MM ****/

      if (PrimitivePolynomial[i] != 0) Alpha_to[MM] ^= mask; /* bit-wise EXOR */
      mask <<= 1;                                        /* single left-shift */
    }
    Index_of[Alpha_to[MM]] = MM;

  /*
   * Have obtained poly-repr of @^MM. Poly-repr of @^(i+1) is given by
   * poly-repr of @^i shifted left one-bit and accounting for any @^MM
   * term that may occur when poly-repr of @^i is shifted.
   */

    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] = A0;
    Alpha_to[NN] = 0;
  }

/*******************************************************************************
*                                                                              *
* Obtain the generator polynomial of the TT-error correcting, length           *
* NN=(2**MM -1) Reed Solomon code from the product of (X+@**(B0+i)), i = 0,    *
* ... ,(2*TT-1)                                                                *
*                                                                              *
* Examples:                                                                    *
*                                                                              *
* If B0 = 1, TT = 1. deg(g(x)) = 2*TT = 2.                                     *
* g(x) = (x+@) (x+@**2)                                                        *
*                                                                              *
* If B0 = 0, TT = 2. deg(g(x)) = 2*TT = 4.                                     *
* g(x) = (x+1) (x+@) (x+@**2) (x+@**3)                                         *
*                                                                              *
*******************************************************************************/

  void generatePolynomial (void) {
    register int i, j;

    GeneratorPolynomial[0] = Alpha_to[B0];
    GeneratorPolynomial[1] = 1;                 /* g(x) = (X+@**B0) initially */
    for (i = 2; i <= NN - KK; i++) {
      GeneratorPolynomial[i] = 1;
/*
 * Below multiply
 *   (GeneratorPolynomial[0]+GeneratorPolynomial[1]*x + ... +GeneratorPolynomial[i]x^i)
 * by (@**(B0+i-1) + x)
 */
      for (j = i - 1; j > 0; j--) {
        if (GeneratorPolynomial[j] != 0)
          GeneratorPolynomial[j] = GeneratorPolynomial[j - 1] ^ Alpha_to[modnn((Index_of[GeneratorPolynomial[j]]) + B0 + i - 1)];
        else
          GeneratorPolynomial[j] = GeneratorPolynomial[j - 1];
      };

    /**** GeneratorPolynomial[0] can never be zero ****/

      GeneratorPolynomial[0] = Alpha_to[modnn((Index_of[GeneratorPolynomial[0]]) + B0 + i - 1)];
    };

  /**** convert GeneratorPolynomial[] to index form for quicker encoding ****/

    for (i = 0; i <= NN - KK; i++) {
      GeneratorPolynomial[i] = Index_of[GeneratorPolynomial[i]];
    };
  }

/*******************************************************************************
*                                                                              *
* take the string of symbols in data[i], i=0..(k-1) and encode                 *
* systematically to produce NN-KK parity symbols in bb[0]..bb[NN-KK-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 GeneratorPolynomial[], which was generated above. Codeword is    *
* c(X) = data(X)*X**(NN-KK)+ b(X)                                              *
*                                                                              *
*******************************************************************************/

  int encode (int data[], int bb[]){
    register int i, j;
    int feedback;

    CLEAR(bb,NN-KK);
    for (i = KK - 1; i >= 0; i--) {
      feedback = Index_of[data[i] ^ bb[NN - KK - 1]];
      if (feedback != A0) {                      /* feedback term is non-zero */
        for (j = NN - KK - 1; j > 0; j--)
          if (GeneratorPolynomial[j] != A0)
            bb[j] = bb[j - 1] ^ Alpha_to[modnn(GeneratorPolynomial[j] + feedback)];
          else
            bb[j] = bb[j - 1];
        bb[0] = Alpha_to[modnn(GeneratorPolynomial[0] + feedback)];
      } else {  /* feedback term is zero, encoder becomes single-byte shifter */
        for (j = NN - KK - 1; j > 0; j--){
          bb[j] = bb[j - 1];
        };
        bb[0] = 0;
      };
    };
    return 0;
  }

/*******************************************************************************
*                                                                              *
* Performs ERRORS+ERASURES decoding of RS codes. If decoding is successful,    *
* writes the codeword into data[] itself. Otherwise data[] is unaltered.       *
*                                                                              *
* Return number of symbols corrected, or -1 if codeword is illegal             *
* or uncorrectable.                                                            *
*                                                                              *
* First "no_eras" erasures are declared by the calling program. Then, the      *
* maximum # of errors correctable is t_after_eras = floor((NN-KK-no_eras)/2).  *
* If the number of channel errors is not greater than "t_after_eras" the       *
* transmitted codeword will be recovered. Details of algorithm can be found    *
* in R. Blahut's "Theory ... of Error-Correcting Codes".                       *
*                                                                              *
*******************************************************************************/

  int decode (int data[], int eras_pos[], int no_eras) {
    int deg_lambda, el, deg_omega;
    int i, j, r;
    int u,q,tmp,num1,num2,den,discr_r;
    int recd[1 << 16];
    int lambda[1 << 16], s[1 << 16];  /* Err+Eras locator & syndrome poly */
    int b[1 << 16], t[1 << 16], omega[1 << 16];
    int root[1 << 16], reg[1 << 16], loc[1 << 16];
    int syn_error, count;

  /**** data[] is in polynomial form, copy and convert to index form ****/

    for (i = NN-1; i >= 0; i--){
      recd[i] = Index_of[data[i]];
    };
/*
 * first form the syndromes; i.e., evaluate recd(x) at roots of g(x)
 * namely @**(B0+i), i = 0, ... ,(NN-KK-1)
 */
    syn_error = 0;
    for (i = 1; i <= NN-KK; i++) {
      tmp = 0;
      for (j = 0; j < NN; j++) {
        if (recd[j] != A0) {                         /* recd[j] in index form */
          tmp ^= Alpha_to[modnn(recd[j] + (B0+i-1)*j)];
        };
      };
      syn_error |= tmp;             /* set flag if non-zero syndrome => error */
      s[i] = Index_of[tmp];                  /* store syndrome in index form  */
    };

    if (!syn_error) {
/*
 * if syndrome is zero, data[] is a codeword and there are no
 * errors to correct. So return data[] unmodified
 */
      return 0;
    };

    CLEAR(&lambda[1],NN-KK);
    lambda[0] = 1;
    if (no_eras > 0) {

    /**** init lambda to be the erasure locator polynomial ****/

      lambda[1] = Alpha_to[eras_pos[0]];
      for (i = 1; i < no_eras; i++) {
        u = 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<NN-KK+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 <= NN-KK) {                            /* 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] != A0)) {
          discr_r ^= Alpha_to[modnn(Index_of[lambda[i]] + s[r - i])];
        };
      };

      discr_r = Index_of[discr_r];                              /* Index form */
      if (discr_r == A0) {
        COPYDOWN(&b[1],b,NN-KK);            /* 2 lines below: B(x) <-- x*B(x) */
        b[0] = A0;
      } else {
        t[0] = lambda[0]; /* 7 lines below: T(x) <-- lambda(x)-discr_r*x*b(x) */
        for (i = 0 ; i < NN-KK; 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 <= NN-KK; i++)
            b[i] = (lambda[i] == 0) ? A0 : modnn(Index_of[lambda[i]] - discr_r + NN);
        } else {
          COPYDOWN(&b[1],b,NN-KK);          /* 2 lines below: B(x) <-- x*B(x) */
          b[0] = A0;
        };
        COPY(lambda,t,NN-KK+1);
      };
    };

  /**** Convert lambda to index form and compute deg(lambda(x)) ****/

    deg_lambda = 0;
    for (i=0;i<NN-KK+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
 */
    COPY(&reg[1],&lambda[1],NN-KK);
    count = 0;                                /* number of roots of lambda(x) */
    for (i = 1; i <= NN; i++) {
      q = 1;
      for (j = deg_lambda; j > 0; j--) {
        if (reg[j] != A0) {
          reg[j] = modnn(reg[j] + j);
          q ^= Alpha_to[reg[j]];
        };
      };
      if (!q) {

      /**** store root (index-form) and error location number ****/

        root[count] = i;
        loc[count] = NN - i;
        count++;
      };
    };

    if (deg_lambda != count) {
/*
 * deg(lambda) unequal to number of roots => uncorrectable
 * error detected
 */
      return -1;
    };

/*
 * Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
 * x**(NN-KK)). in index form. Also find deg(omega).
 */
    deg_omega = 0;
    for (i = 0; i < NN-KK;i++){
      tmp = 0;
      j = (deg_lambda < i) ? deg_lambda : i;
      for (;j >= 0; j--){
        if ((s[i + 1 - j] != A0) && (lambda[j] != A0))
          tmp ^= Alpha_to[modnn(s[i + 1 - j] + lambda[j])];
      };
      if (tmp != 0) deg_omega = i;
      omega[i] = Index_of[tmp];
    };
    omega[NN-KK] = A0;

/*
 * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
 * inv(X(l))**(B0-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] * (B0 - 1) + NN)];
      den = 0;

    /**** lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] ****/

      for (i = min(deg_lambda,NN-KK-1) & ~1; i >= 0; i -=2) {
        if (lambda[i+1] != A0)
          den ^= Alpha_to[modnn(lambda[i+1] + i * root[j])];
      };

      if (den == 0) {
        return -1;
      };

    /**** Apply error to data ****/

      if (num1 != 0) {
        data[loc[j]] ^= Alpha_to[modnn(Index_of[num1] + Index_of[num2] + NN - Index_of[den])];
      };
    };
    return count;
  }