home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / gnu / ecc-1.2.1 / rslib.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-01-19  |  8.2 KB  |  334 lines
  1. /*
  2.     ecc Version 1.2  by Paul Flaherty (paulf@stanford.edu)
  3.     Copyright (C) 1993 Free Software Foundation, Inc.
  4.  
  5.     This program is free software; you can redistribute it and/or modify
  6.     it under the terms of the GNU General Public License as published by
  7.     the Free Software Foundation; either version 2, or (at your option)
  8.     any later version.
  9.  
  10.     This program is distributed in the hope that it will be useful,
  11.     but WITHOUT ANY WARRANTY; without even the implied warranty of
  12.     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13.     GNU General Public License for more details.
  14.  
  15.     You should have received a copy of the GNU General Public License
  16.     along with this program; if not, write to the Free Software
  17.     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18. */
  19.  
  20.  
  21. /* rslib.c
  22.     Library of Reed - Solomon Routines
  23.  
  24.     This file contains the actual routines to implement a Reed - Solomon
  25.     (255,249,7) code.  The encoder uses a feedback shift register
  26.     generator, which systematically encodes 249 bytes into a 255 byte
  27.     block.  The decoder is a classic Peterson algorithm.
  28. */
  29.  
  30. #include "ecc.h"
  31.  
  32.  
  33. /* Reed - Solomon Encoder.  The Encoder uses a shift register algorithm,
  34.    as detailed in _Applied Modern Algebra_ by Dornhoff and Hohn (p.446).
  35.    Note that the message is reversed in the code array; this was done to
  36.    allow for (emergency) recovery of the message directly from the
  37.    data stream.
  38. */
  39.  
  40. rsencode (m, c)
  41.  
  42.      unsigned char m[249], c[255];
  43.  
  44. {
  45.  
  46.   unsigned char r[6], rtmp;
  47.  
  48.   int i, j;
  49.  
  50.   for (i = 0; i < 6; i++)
  51.     r[i] = 0;
  52.  
  53.   for (i = 0; i < 249; i++)
  54.     {
  55.       c[254 - i] = m[i];
  56.       rtmp = gfadd (m[i], r[5]);
  57.       for (j = 5; j > 0; j--)
  58.     {
  59.       r[j] = gfadd (gfmul (rtmp, g[j]), r[j - 1]);
  60.     }
  61.       r[0] = gfmul (rtmp, g[0]);
  62.     }
  63.   for (i = 0; i < 6; i++)
  64.     {
  65.       c[i] = r[i];
  66.     }
  67. }
  68.  
  69.  
  70. /* Polynomial Evaluator, used to determine the Syndrome Vector.  This is
  71.    relatively straightforward, and there are faster algorithms.
  72. */
  73.  
  74. unsigned char
  75. evalpoly (p, x)
  76.  
  77.      unsigned char p[255], x;
  78.  
  79. {
  80.  
  81.   unsigned char y;
  82.   int i;
  83.   y = 0;
  84.  
  85.   for (i = 0; i < 255; i++)
  86.     {
  87.       y = gfadd (y, gfmul (p[i], gfexp (x, i)));
  88.     }
  89.   return (y);
  90.  
  91. }
  92.  
  93.  
  94. /* Determine the Syndrome Vector.  Note that in s[0] we return the OR of
  95.    all of the syndromes; this allows for an easy check for the no - error
  96.    condition.
  97. */
  98.  
  99. syndrome (c, s)
  100.  
  101.      unsigned char c[255], s[7];
  102. {
  103.  
  104.   extern unsigned char e2v[256];
  105.   int i;
  106.   s[0] = 0;
  107.   for (i = 1; i < 7; i++)
  108.     {
  109.       s[i] = evalpoly (c, e2v[i]);
  110.       s[0] = s[0] | s[i];
  111.     }
  112.  
  113. }
  114.  
  115.  
  116. /* Determine the number of errors in a block.  Since we have to find the
  117.    determinant of the S[] matrix in order to determine singularity, we
  118.    also return the determinant to be used by the Cramer's Rule correction
  119.    algorithm.
  120. */
  121.  
  122. errnum (s, det, errs)
  123.  
  124.      unsigned char s[7], *det;
  125.      int *errs;
  126.  
  127. {
  128.  
  129.  
  130.   *det = gfmul (s[2], gfmul (s[4], s[6]));
  131.   *det = gfadd (*det, gfmul (s[2], gfmul (s[5], s[5])));
  132.   *det = gfadd (*det, gfmul (s[6], gfmul (s[3], s[3])));
  133.   *det = gfadd (*det, gfmul (s[4], gfmul (s[4], s[4])));
  134.   *errs = 3;
  135.  
  136.   if (*det != 0)
  137.     return;
  138.  
  139.   *det = gfadd (gfmul (s[2], s[4]), gfexp (s[3], 2));
  140.   *errs = 2;
  141.   if (*det != 0)
  142.     return;
  143.  
  144.   *det = s[1];
  145.   *errs = 1;
  146.   if (*det != 0)
  147.     return;
  148.  
  149.   *errs = 4;
  150.  
  151. }
  152.  
  153.  
  154. /* Full impementation of the three error correcting Peterson decoder.  For
  155.    t<6, it is faster than Massey - Berlekamp.  It is also somewhat more
  156.    intuitive.
  157. */
  158.  
  159. rsdecode (code, mesg, errcode)
  160.  
  161.      unsigned char code[255], mesg[249];
  162.      int *errcode;
  163. {
  164.  
  165.   extern unsigned char v2e[256];
  166.   unsigned char syn[7], deter, z[4], e0, e1, e2, n0, n1, n2, w0, w1, w2,
  167.     x0, x[3];
  168.   int i, sols;
  169.  
  170.   *errcode = 0;
  171.  
  172.   /* First, get the message out of the code, so that even if we can't correct
  173.        it, we return an estimate.
  174.     */
  175.  
  176.   for (i = 0; i < 249; i++)
  177.     mesg[i] = code[254 - i];
  178.  
  179.   syndrome (code, syn);
  180.  
  181.   if (syn[0] == 0)
  182.     return;
  183.  
  184.   /* We now know we have at least one error.  If there are no errors detected,
  185.        we assume that something funny is going on, and so return with errcode 4,
  186.        else pass the number of errors back via errcode.
  187.     */
  188.  
  189.   errnum (syn, &deter, errcode);
  190.  
  191.   if (*errcode == 4)
  192.     return;
  193.  
  194.   /* Having obtained the syndrome, the number of errors, and the determinant,
  195.        we now proceed to correct the block.  If we do not find exactly the
  196.        number of solutions equal to the number of errors, we have exceeded our
  197.        error capacity, and return with the block uncorrected, and errcode 4.
  198.     */
  199.  
  200.   switch (*errcode)
  201.     {
  202.  
  203.     case 1:
  204.  
  205.       x0 = gfmul (syn[2], gfinv (syn[1]));
  206.       w0 = gfmul (gfexp (syn[1], 2), gfinv (syn[2]));
  207.       if (v2e[x0] > 5)
  208.     mesg[254 - v2e[x0]] = gfadd (mesg[254 - v2e[x0]], w0);
  209.       return;
  210.  
  211.     case 2:
  212.  
  213.       z[0] = gfmul (gfadd (gfmul (syn[1], syn[3]), gfexp (syn[2], 2)), gfinv (deter));
  214.       z[1] = gfmul (gfadd (gfmul (syn[2], syn[3]), gfmul (syn[1], syn[4])), gfinv (deter));
  215.       z[2] = 1;
  216.       z[3] = 0;
  217.  
  218.       polysolve (z, x, &sols);
  219.       if (sols != 2)
  220.     {
  221.       *errcode = 4;
  222.       return;
  223.     }
  224.  
  225.       w0 = gfmul (z[0], syn[1]);
  226.       w1 = gfadd (gfmul (z[0], syn[2]), gfmul (z[1], syn[1]));
  227.       n0 = 254 - v2e[gfinv (x[0])];
  228.       n1 = 254 - v2e[gfinv (x[1])];
  229.       e0 = gfmul (gfadd (w0, gfmul (w1, x[0])), gfinv (z[1]));
  230.       e1 = gfmul (gfadd (w0, gfmul (w1, x[1])), gfinv (z[1]));
  231.  
  232.       if (n0 < 249)
  233.     mesg[n0] = gfadd (mesg[n0], e0);
  234.       if (n1 < 249)
  235.     mesg[n1] = gfadd (mesg[n1], e1);
  236.       return;
  237.  
  238.     case 3:
  239.  
  240.       z[3] = 1;
  241.       z[2] = gfmul (syn[1], gfmul (syn[4], syn[6]));
  242.       z[2] = gfadd (z[2], gfmul (syn[1], gfmul (syn[5], syn[5])));
  243.       z[2] = gfadd (z[2], gfmul (syn[5], gfmul (syn[3], syn[3])));
  244.       z[2] = gfadd (z[2], gfmul (syn[3], gfmul (syn[4], syn[4])));
  245.       z[2] = gfadd (z[2], gfmul (syn[2], gfmul (syn[5], syn[4])));
  246.       z[2] = gfadd (z[2], gfmul (syn[2], gfmul (syn[3], syn[6])));
  247.       z[2] = gfmul (z[2], gfinv (deter));
  248.  
  249.       z[1] = gfmul (syn[1], gfmul (syn[3], syn[6]));
  250.       z[1] = gfadd (z[1], gfmul (syn[1], gfmul (syn[5], syn[4])));
  251.       z[1] = gfadd (z[1], gfmul (syn[4], gfmul (syn[3], syn[3])));
  252.       z[1] = gfadd (z[1], gfmul (syn[2], gfmul (syn[4], syn[4])));
  253.       z[1] = gfadd (z[1], gfmul (syn[2], gfmul (syn[3], syn[5])));
  254.       z[1] = gfadd (z[1], gfmul (syn[2], gfmul (syn[2], syn[6])));
  255.       z[1] = gfmul (z[1], gfinv (deter));
  256.  
  257.       z[0] = gfmul (syn[2], gfmul (syn[3], syn[4]));
  258.       z[0] = gfadd (z[0], gfmul (syn[3], gfmul (syn[2], syn[4])));
  259.       z[0] = gfadd (z[0], gfmul (syn[3], gfmul (syn[5], syn[1])));
  260.       z[0] = gfadd (z[0], gfmul (syn[4], gfmul (syn[4], syn[1])));
  261.       z[0] = gfadd (z[0], gfmul (syn[3], gfmul (syn[3], syn[3])));
  262.       z[0] = gfadd (z[0], gfmul (syn[2], gfmul (syn[2], syn[5])));
  263.       z[0] = gfmul (z[0], gfinv (deter));
  264.  
  265.       polysolve (z, x, &sols);
  266.       if (sols != 3)
  267.     {
  268.       *errcode = 4;
  269.       return;
  270.     }
  271.  
  272.       w0 = gfmul (z[0], syn[1]);
  273.       w1 = gfadd (gfmul (z[0], syn[2]), gfmul (z[1], syn[1]));
  274.       w2 = gfadd (gfmul (z[0], syn[3]), gfadd (gfmul (z[1], syn[2]), gfmul (z[2], syn[1])));
  275.  
  276.       n0 = 254 - v2e[gfinv (x[0])];
  277.       n1 = 254 - v2e[gfinv (x[1])];
  278.       n2 = 254 - v2e[gfinv (x[2])];
  279.  
  280.       e0 = gfadd (w0, gfadd (gfmul (w1, x[0]), gfmul (w2, gfexp (x[0], 2))));
  281.       e0 = gfmul (e0, gfinv (gfadd (z[1], gfexp (x[0], 2))));
  282.       e1 = gfadd (w0, gfadd (gfmul (w1, x[1]), gfmul (w2, gfexp (x[1], 2))));
  283.       e1 = gfmul (e1, gfinv (gfadd (z[1], gfexp (x[1], 2))));
  284.       e2 = gfadd (w0, gfadd (gfmul (w1, x[2]), gfmul (w2, gfexp (x[2], 2))));
  285.       e2 = gfmul (e2, gfinv (gfadd (z[1], gfexp (x[2], 2))));
  286.  
  287.       if (n0 < 249)
  288.     mesg[n0] = gfadd (mesg[n0], e0);
  289.       if (n1 < 249)
  290.     mesg[n1] = gfadd (mesg[n1], e1);
  291.       if (n2 < 249)
  292.     mesg[n2] = gfadd (mesg[n2], e2);
  293.  
  294.       return;
  295.  
  296.     default:
  297.  
  298.       *errcode = 4;
  299.       return;
  300.  
  301.     }
  302.  
  303. }
  304.  
  305.  
  306. /* Polynomial Solver.  Simple exhaustive search, as solving polynomials is
  307.    generally NP - Complete anyway.
  308. */
  309.  
  310. polysolve (polynom, roots, numsol)
  311.  
  312.      unsigned char polynom[4], roots[3];
  313.      int *numsol;
  314.  
  315. {
  316.   extern unsigned char e2v[256];
  317.   int i, j;
  318.   unsigned char y;
  319.  
  320.   *numsol = 0;
  321.  
  322.   for (i = 0; i < 255; i++)
  323.     {
  324.       y = 0;
  325.       for (j = 0; j < 4; j++)
  326.     y = gfadd (y, gfmul (polynom[j], gfexp (e2v[i], j)));
  327.       if (y == 0)
  328.     {
  329.       roots[*numsol] = e2v[i];
  330.       *numsol = *numsol + 1;
  331.     }
  332.     }
  333. }
  334.