"Fossies" - the Fresh Open Source Software Archive

Member "dvdisaster-0.79.5/galois.c" (3 Jan 2015, 5752 Bytes) of package /linux/misc/dvdisaster-0.79.5.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "galois.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes reports: 0.79.3_vs_0.79.5 or 0.72.6_vs_0.79.5.

    1 /*  dvdisaster: Additional error correction for optical media.
    2  *  Copyright (C) 2004-2015 Carsten Gnoerlich.
    3  *
    4  *  The Reed-Solomon error correction draws a lot of inspiration - and even code -
    5  *  from Phil Karn's excellent Reed-Solomon library: http://www.ka9q.net/code/fec/
    6  *
    7  *  Email: carsten@dvdisaster.org  -or-  cgnoerlich@fsfe.org
    8  *  Project homepage: http://www.dvdisaster.org
    9  *
   10  *  This file is part of dvdisaster.
   11  *
   12  *  dvdisaster is free software: you can redistribute it and/or modify
   13  *  it under the terms of the GNU General Public License as published by
   14  *  the Free Software Foundation, either version 3 of the License, or
   15  *  (at your option) any later version.
   16  *
   17  *  dvdisaster is distributed in the hope that it will be useful,
   18  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   19  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   20  *  GNU General Public License for more details.
   21  *
   22  *  You should have received a copy of the GNU General Public License
   23  *  along with dvdisaster. If not, see <http://www.gnu.org/licenses/>.
   24  */
   25 
   26 #include "dvdisaster.h"
   27 
   28 #include "galois-inlines.h"
   29 
   30 /***
   31  *** Galois field arithmetic.
   32  *** 
   33  * Calculations are done over the extension field GF(2**n).
   34  * Be careful not to overgeneralize these arithmetics;
   35  * they only work for the case of GF(p**n) with p being prime.
   36  */
   37 
   38 /* Initialize the Galois field tables */
   39 
   40 
   41 GaloisTables* CreateGaloisTables(gint32 gf_generator)
   42 {  GaloisTables *gt = g_malloc0(sizeof(GaloisTables));
   43    gint32 b,log;
   44 
   45    /* Allocate the tables.
   46       The encoder uses a special version of alpha_to which has the mod_fieldmax()
   47       folded into the table. */
   48 
   49    gt->gfGenerator = gf_generator;
   50 
   51    gt->indexOf     = g_malloc(GF_FIELDSIZE * sizeof(gint32));
   52    gt->alphaTo     = g_malloc(GF_FIELDSIZE * sizeof(gint32));
   53    gt->encAlphaTo  = g_malloc(2*GF_FIELDSIZE * sizeof(gint32));
   54    
   55    /* create the log/ilog values */
   56 
   57    for(b=1, log=0; log<GF_FIELDMAX; log++)
   58    {  gt->indexOf[b]   = log;
   59       gt->alphaTo[log] = b;
   60       b = b << 1;
   61       if(b & GF_FIELDSIZE)
   62     b = b ^ gf_generator;
   63    }
   64 
   65    if(b!=1) Stop("Failed to create the Galois field log tables!\n");
   66 
   67    /* we're even closed using infinity (makes things easier) */
   68 
   69    gt->indexOf[0] = GF_ALPHA0;    /* log(0) = inf */
   70    gt->alphaTo[GF_ALPHA0] = 0;   /* and the other way around */
   71 
   72    for(b=0; b<2*GF_FIELDSIZE; b++)
   73      gt->encAlphaTo[b] = gt->alphaTo[mod_fieldmax(b)];
   74 
   75    return gt;
   76 }
   77 
   78 void FreeGaloisTables(GaloisTables *gt)
   79 {
   80   if(gt->indexOf)     g_free(gt->indexOf);
   81   if(gt->alphaTo)     g_free(gt->alphaTo);
   82   if(gt->encAlphaTo) g_free(gt->encAlphaTo);
   83 
   84   g_free(gt);
   85 }
   86 
   87 /***
   88  *** Create the Reed-Solomon generator polynomial
   89  *** and some auxiliary data structures.
   90  */
   91 
   92 ReedSolomonTables *CreateReedSolomonTables(GaloisTables *gt,
   93                        gint32 first_consecutive_root,
   94                        gint32 prim_elem,
   95                        int nroots_in)
   96 {  ReedSolomonTables *rt = g_malloc0(sizeof(ReedSolomonTables));
   97    int lut_size, feedback;
   98    gint32 i,j,root;
   99    guint8 *lut;
  100 
  101    rt->gfTables = gt;
  102    rt->fcr      = first_consecutive_root;
  103    rt->primElem = prim_elem;
  104    rt->nroots   = nroots_in;
  105    rt->ndata    = GF_FIELDMAX - rt->nroots;
  106 
  107    rt->gpoly    = g_malloc((rt->nroots+1) * sizeof(gint32));
  108 
  109    /* Create the RS code generator polynomial */
  110 
  111    rt->gpoly[0] = 1;
  112 
  113    for(i=0, root=first_consecutive_root*prim_elem; i<rt->nroots; i++, root+=prim_elem)
  114    {  rt->gpoly[i+1] = 1;
  115 
  116      /* Multiply gpoly  by  alpha**(root+x) */
  117 
  118      for(j=i; j>0; j--)
  119      {
  120        if(rt->gpoly[j] != 0)
  121          rt->gpoly[j] = rt->gpoly[j-1] ^ gt->alphaTo[mod_fieldmax(gt->indexOf[rt->gpoly[j]] + root)]; 
  122        else
  123      rt->gpoly[j] = rt->gpoly[j-1];
  124      }
  125 
  126      rt->gpoly[0] = gt->alphaTo[mod_fieldmax(gt->indexOf[rt->gpoly[0]] + root)];
  127    }
  128 
  129    /* Store the polynomials index for faster encoding */ 
  130 
  131    for(i=0; i<=rt->nroots; i++)
  132      rt->gpoly[i] = gt->indexOf[rt->gpoly[i]];
  133 
  134 #if 0
  135    /* for the precalculated unrolled loops only */
  136 
  137    for(i=gt->nroots-1; i>0; i--)
  138      PrintCLI(
  139         "                  par_idx[((++spk)&%d)] ^= enc_alpha_to[feedback + %3d];\n",
  140         nroots-1,gt->gpoly[i]);
  141 
  142    PrintCLI("                  par_idx[sp] = enc_alpha_to[feedback + %3d];\n",
  143       gt->gpoly[0]);
  144 #endif
  145 
  146   /* 
  147    * Initialize the shift pointer so that we will come out at shiftPtr==0
  148    * respectively (ndata+sp) mod nroots = 0 after working in all ndata layers.
  149    */
  150 
  151    rt->shiftInit = rt->nroots - rt->ndata % rt->nroots;
  152    if(rt->shiftInit == rt->nroots)
  153      rt->shiftInit = 0;
  154 
  155    /*
  156     * Initialize lookup tables for both encoder types.
  157     * The 32bit portable encoder will shift them to word boundaries,
  158     * while the SSE2 encoder does direct unaligned reads.
  159     */
  160 
  161    lut_size = (rt->nroots+15)&~15;
  162    lut_size += 16;
  163    for(i=0; i<GF_FIELDSIZE; i++)
  164       rt->bLut[i] = g_malloc0(2*lut_size);
  165 
  166    for(feedback=0; feedback<256; feedback++)
  167    {  gint32 *gpoly        = rt->gpoly + rt->nroots;
  168       gint32 *enc_alpha_to = gt->encAlphaTo;
  169       int nroots = rt->nroots;
  170 
  171       for(i=0; i<nroots; i++)
  172       {  guint8 value = (guint8)enc_alpha_to[feedback + *--gpoly];
  173      rt->bLut[feedback][i] = rt->bLut[feedback][nroots+i] = value; 
  174       }
  175    }
  176 
  177    /*
  178     * Prepare lookup table for syndrome calculation.
  179     */
  180 
  181    lut = rt->synLut = g_malloc(rt->nroots * GF_FIELDSIZE * sizeof(int));
  182    for(i=0; i<rt->nroots; i++)
  183      for(j=0; j<GF_FIELDSIZE; j++)
  184        *lut++ = gt->alphaTo[mod_fieldmax(gt->indexOf[j] + (rt->fcr+i)*rt->primElem)];
  185 
  186    return rt;
  187 }
  188 
  189 void FreeReedSolomonTables(ReedSolomonTables *rt)
  190 { int i;
  191 
  192   if(rt->gpoly)        g_free(rt->gpoly);
  193 
  194   for(i=0; i<GF_FIELDSIZE; i++)
  195   {  g_free(rt->bLut[i]);
  196   }
  197   g_free(rt->synLut);
  198 
  199   g_free(rt);
  200 }