bch15_5.c

00001 /*Copyright (C) 2008-2009  Timothy B. Terriberry (tterribe@xiph.org)
00002   You can redistribute this library and/or modify it under the terms of the
00003    GNU Lesser General Public License as published by the Free Software
00004    Foundation; either version 2.1 of the License, or (at your option) any later
00005    version.*/
00006 #include "bch15_5.h"
00007 
00008 /*A cycle in GF(2**4) generated by alpha=(x**4+x+1).
00009   It is extended an extra 16 entries to avoid some expensive mod operations.*/
00010 static const unsigned char gf16_exp[31]={
00011   1,2,4,8,3,6,12,11,5,10,7,14,15,13,9,1,2,4,8,3,6,12,11,5,10,7,14,15,13,9,1
00012 };
00013 
00014 /*The location of each integer 1...16 in the cycle.*/
00015 static const signed char gf16_log[16]={
00016   -1,0,1,4,2,8,5,10,3,14,9,7,6,13,11,12
00017 };
00018 
00019 /*Multiplication in GF(2**4) using logarithms.*/
00020 static unsigned gf16_mul(unsigned _a,unsigned _b){
00021   return _a==0||_b==0?0:gf16_exp[gf16_log[_a]+gf16_log[_b]];
00022 }
00023 
00024 /*Division in GF(2**4) using logarithms.
00025   The result when dividing by zero is undefined.*/
00026 static unsigned gf16_div(unsigned _a,unsigned _b){
00027   return _a==0?0:gf16_exp[gf16_log[_a]+15-gf16_log[_b]];
00028 }
00029 
00030 /*Multiplication in GF(2**4) when the second argument is known to be non-zero
00031    (proven by representing it by its logarithm).*/
00032 static unsigned gf16_hmul(unsigned _a,unsigned _logb){
00033   return _a==0?0:gf16_exp[gf16_log[_a]+_logb];
00034 }
00035 
00036 /*The syndrome normally has five values, S_1 ... S_5.
00037   We only calculate and store the odd ones in _s, since S_2=S_1**2 and
00038    S_4=S_2**2.
00039   Returns zero iff all the syndrome values are zero.*/
00040 static int bch15_5_calc_syndrome(unsigned _s[3],unsigned _y){
00041   unsigned p;
00042   int      i;
00043   int      j;
00044   p=0;
00045   for(i=0;i<15;i++)if(_y&1<<i)p^=gf16_exp[i];
00046   _s[0]=p;
00047   p=0;
00048   for(i=0;i<3;i++)for(j=0;j<5;j++)if(_y&1<<5*i+j)p^=gf16_exp[j*3];
00049   _s[1]=p;
00050   p=0;
00051   for(i=0;i<5;i++)for(j=0;j<3;j++)if(_y&1<<3*i+j)p^=gf16_exp[j*5];
00052   _s[2]=p;
00053   return _s[0]!=0||_s[1]!=0||_s[2]!=0;
00054 }
00055 
00056 /*Compute the coefficients of the error-locator polynomial.
00057   Returns the number of errors (the degree of the polynomial).*/
00058 static int bch15_5_calc_omega(unsigned _o[3],unsigned _s[3]){
00059   unsigned s02;
00060   unsigned tt;
00061   unsigned dd;
00062   int      d;
00063   _o[0]=_s[0];
00064   s02=gf16_mul(_s[0],_s[0]);
00065   dd=_s[1]^gf16_mul(_s[0],s02);
00066   tt=_s[2]^gf16_mul(s02,_s[1]);
00067   _o[1]=dd?gf16_div(tt,dd):0;
00068   _o[2]=dd^gf16_mul(_s[0],_o[1]);
00069   for(d=3;d>0&&!_o[d-1];d--);
00070   return d;
00071 }
00072 
00073 /*Find the roots of the error polynomial.
00074   Returns the number of roots found, or a negative value if the polynomial did
00075    not have enough roots, indicating a decoding error.*/
00076 static int bch15_5_calc_epos(unsigned _epos[3],unsigned _s[3]){
00077   unsigned o[3];
00078   int      nerrors;
00079   int      d;
00080   int      i;
00081   d=bch15_5_calc_omega(o,_s);
00082   nerrors=0;
00083   if(d==1)_epos[nerrors++]=gf16_log[o[0]];
00084   else if(d>0){
00085     for(i=0;i<15;i++){
00086       int i2;
00087       i2=gf16_log[gf16_exp[i<<1]];
00088       if(!(gf16_exp[i+i2]^gf16_hmul(o[0],i2)^gf16_hmul(o[1],i)^o[2])){
00089         _epos[nerrors++]=i;
00090       }
00091     }
00092     if(nerrors<d)return -1;
00093   }
00094   return nerrors;
00095 }
00096 
00097 int bch15_5_correct(unsigned *_y){
00098   unsigned s[3];
00099   unsigned epos[3];
00100   unsigned y;
00101   int      nerrors;
00102   int      i;
00103   y=*_y;
00104   if(!bch15_5_calc_syndrome(s,y))return 0;
00105   nerrors=bch15_5_calc_epos(epos,s);
00106   if(nerrors>0){
00107     /*If we had a non-zero syndrome value, we should always find at least one
00108        error location, or we've got a decoding error.*/
00109     for(i=0;i<nerrors;i++)y^=1<<epos[i];
00110     /*If there were too many errors, we may not find enough roots to reduce the
00111        syndrome to zero.
00112       We could recompute it to check, but it's much faster just to check that
00113        we have a valid codeword.*/
00114     if(bch15_5_encode(y>>10)==y){
00115       /*Decoding succeeded.*/
00116       *_y=y;
00117       return nerrors;
00118     }
00119   }
00120   /*Decoding failed due to too many bit errors.*/
00121   return -1;
00122 }
00123 
00124 unsigned bch15_5_encode(unsigned _x){
00125   return (-(_x&1)&0x0537)^(-(_x>>1&1)&0x0A6E)^(-(_x>>2&1)&0x11EB)^
00126    (-(_x>>3&1)&0x23D6)^(-(_x>>4&1)&0x429B);
00127 }
00128 
00129 #if 0
00130 #include <stdio.h>
00131 
00132 static unsigned codes[32];
00133 
00134 static int hamming(int _a,int _b){
00135   int d;
00136   int n;
00137   d=_a^_b;
00138   for(n=0;d;n++)d&=d-1;
00139   return n;
00140 }
00141 
00142 static int closest(int _y){
00143   int min_i;
00144   int min_d;
00145   int i;
00146   int d;
00147   min_i=0;
00148   min_d=hamming(_y,codes[0]);
00149   for(i=1;i<32;i++){
00150     d=hamming(_y,codes[i]);
00151     if(d<min_d){
00152       min_d=d;
00153       min_i=i;
00154     }
00155   }
00156   return codes[min_i];
00157 }
00158 
00159 int main(void){
00160   int i;
00161   /*Print a list of the valid (uncorrupt) codewords.*/
00162   for(i=0;i<32;i++)codes[i]=bch15_5_encode(i);
00163   for(i=0;i<32;i++)printf("0x%04X%s",codes[i],i+1<32?"  ":"\n");
00164   /*Try to decode all receivable (possibly corrupt) codewords.*/
00165   for(i=0;i<0x8000;i++){
00166     unsigned y;
00167     unsigned z;
00168     int      nerrors;
00169     int      j;
00170     y=i;
00171     nerrors=bch15_5_correct(&y);
00172     z=closest(i);
00173     if(nerrors<0){
00174       printf("0x%04X->Failed\n",i);
00175       if(hamming(i,z)<=3)printf("Error: 0x%04X should map to 0x%04X\n",i,z);
00176     }
00177     else{
00178       printf("0x%04X->0x%04X\n",i,y);
00179       if(z!=y)printf("Error: 0x%04X should map to 0x%04X\n",i,z);
00180     }
00181   }
00182   return 0;
00183 }
00184 #endif
00185