185 lines
4.9 KiB
C
185 lines
4.9 KiB
C
/*Copyright (C) 2008-2009 Timothy B. Terriberry (tterribe@xiph.org)
|
|
You can redistribute this library and/or modify it under the terms of the
|
|
GNU Lesser General Public License as published by the Free Software
|
|
Foundation; either version 2.1 of the License, or (at your option) any later
|
|
version.*/
|
|
#include "bch15_5.h"
|
|
|
|
/*A cycle in GF(2**4) generated by alpha=(x**4+x+1).
|
|
It is extended an extra 16 entries to avoid some expensive mod operations.*/
|
|
static const unsigned char gf16_exp[31]={
|
|
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
|
|
};
|
|
|
|
/*The location of each integer 1...16 in the cycle.*/
|
|
static const signed char gf16_log[16]={
|
|
-1,0,1,4,2,8,5,10,3,14,9,7,6,13,11,12
|
|
};
|
|
|
|
/*Multiplication in GF(2**4) using logarithms.*/
|
|
static unsigned gf16_mul(unsigned _a,unsigned _b){
|
|
return _a==0||_b==0?0:gf16_exp[gf16_log[_a]+gf16_log[_b]];
|
|
}
|
|
|
|
/*Division in GF(2**4) using logarithms.
|
|
The result when dividing by zero is undefined.*/
|
|
static unsigned gf16_div(unsigned _a,unsigned _b){
|
|
return _a==0?0:gf16_exp[gf16_log[_a]+15-gf16_log[_b]];
|
|
}
|
|
|
|
/*Multiplication in GF(2**4) when the second argument is known to be non-zero
|
|
(proven by representing it by its logarithm).*/
|
|
static unsigned gf16_hmul(unsigned _a,unsigned _logb){
|
|
return _a==0?0:gf16_exp[gf16_log[_a]+_logb];
|
|
}
|
|
|
|
/*The syndrome normally has five values, S_1 ... S_5.
|
|
We only calculate and store the odd ones in _s, since S_2=S_1**2 and
|
|
S_4=S_2**2.
|
|
Returns zero iff all the syndrome values are zero.*/
|
|
static int bch15_5_calc_syndrome(unsigned _s[3],unsigned _y){
|
|
unsigned p;
|
|
int i;
|
|
int j;
|
|
p=0;
|
|
for(i=0;i<15;i++)if(_y&1<<i)p^=gf16_exp[i];
|
|
_s[0]=p;
|
|
p=0;
|
|
for(i=0;i<3;i++)for(j=0;j<5;j++)if(_y&1<<(5*i+j))p^=gf16_exp[j*3];
|
|
_s[1]=p;
|
|
p=0;
|
|
for(i=0;i<5;i++)for(j=0;j<3;j++)if(_y&1<<(3*i+j))p^=gf16_exp[j*5];
|
|
_s[2]=p;
|
|
return _s[0]!=0||_s[1]!=0||_s[2]!=0;
|
|
}
|
|
|
|
/*Compute the coefficients of the error-locator polynomial.
|
|
Returns the number of errors (the degree of the polynomial).*/
|
|
static int bch15_5_calc_omega(unsigned _o[3],unsigned _s[3]){
|
|
unsigned s02;
|
|
unsigned tt;
|
|
unsigned dd;
|
|
int d;
|
|
_o[0]=_s[0];
|
|
s02=gf16_mul(_s[0],_s[0]);
|
|
dd=_s[1]^gf16_mul(_s[0],s02);
|
|
tt=_s[2]^gf16_mul(s02,_s[1]);
|
|
_o[1]=dd?gf16_div(tt,dd):0;
|
|
_o[2]=dd^gf16_mul(_s[0],_o[1]);
|
|
for(d=3;d>0&&!_o[d-1];d--);
|
|
return d;
|
|
}
|
|
|
|
/*Find the roots of the error polynomial.
|
|
Returns the number of roots found, or a negative value if the polynomial did
|
|
not have enough roots, indicating a decoding error.*/
|
|
static int bch15_5_calc_epos(unsigned _epos[3],unsigned _s[3]){
|
|
unsigned o[3];
|
|
int nerrors;
|
|
int d;
|
|
int i;
|
|
d=bch15_5_calc_omega(o,_s);
|
|
nerrors=0;
|
|
if(d==1)_epos[nerrors++]=gf16_log[o[0]];
|
|
else if(d>0){
|
|
for(i=0;i<15;i++){
|
|
int i2;
|
|
i2=gf16_log[gf16_exp[i<<1]];
|
|
if(!(gf16_exp[i+i2]^gf16_hmul(o[0],i2)^gf16_hmul(o[1],i)^o[2])){
|
|
_epos[nerrors++]=i;
|
|
}
|
|
}
|
|
if(nerrors<d)return -1;
|
|
}
|
|
return nerrors;
|
|
}
|
|
|
|
int bch15_5_correct(unsigned *_y){
|
|
unsigned s[3];
|
|
unsigned epos[3];
|
|
unsigned y;
|
|
int nerrors;
|
|
int i;
|
|
y=*_y;
|
|
if(!bch15_5_calc_syndrome(s,y))return 0;
|
|
nerrors=bch15_5_calc_epos(epos,s);
|
|
if(nerrors>0){
|
|
/*If we had a non-zero syndrome value, we should always find at least one
|
|
error location, or we've got a decoding error.*/
|
|
for(i=0;i<nerrors;i++)y^=1<<epos[i];
|
|
/*If there were too many errors, we may not find enough roots to reduce the
|
|
syndrome to zero.
|
|
We could recompute it to check, but it's much faster just to check that
|
|
we have a valid codeword.*/
|
|
if(bch15_5_encode(y>>10)==y){
|
|
/*Decoding succeeded.*/
|
|
*_y=y;
|
|
return nerrors;
|
|
}
|
|
}
|
|
/*Decoding failed due to too many bit errors.*/
|
|
return -1;
|
|
}
|
|
|
|
unsigned bch15_5_encode(unsigned _x){
|
|
return (-(_x&1)&0x0537)^(-(_x>>1&1)&0x0A6E)^(-(_x>>2&1)&0x11EB)^
|
|
(-(_x>>3&1)&0x23D6)^(-(_x>>4&1)&0x429B);
|
|
}
|
|
|
|
#if 0
|
|
#include <stdio.h>
|
|
|
|
static unsigned codes[32];
|
|
|
|
static int hamming(int _a,int _b){
|
|
int d;
|
|
int n;
|
|
d=_a^_b;
|
|
for(n=0;d;n++)d&=d-1;
|
|
return n;
|
|
}
|
|
|
|
static int closest(int _y){
|
|
int min_i;
|
|
int min_d;
|
|
int i;
|
|
int d;
|
|
min_i=0;
|
|
min_d=hamming(_y,codes[0]);
|
|
for(i=1;i<32;i++){
|
|
d=hamming(_y,codes[i]);
|
|
if(d<min_d){
|
|
min_d=d;
|
|
min_i=i;
|
|
}
|
|
}
|
|
return codes[min_i];
|
|
}
|
|
|
|
int main(void){
|
|
int i;
|
|
/*Print a list of the valid (uncorrupt) codewords.*/
|
|
for(i=0;i<32;i++)codes[i]=bch15_5_encode(i);
|
|
for(i=0;i<32;i++)printf("0x%04X%s",codes[i],i+1<32?" ":"\n");
|
|
/*Try to decode all receivable (possibly corrupt) codewords.*/
|
|
for(i=0;i<0x8000;i++){
|
|
unsigned y;
|
|
unsigned z;
|
|
int nerrors;
|
|
int j;
|
|
y=i;
|
|
nerrors=bch15_5_correct(&y);
|
|
z=closest(i);
|
|
if(nerrors<0){
|
|
printf("0x%04X->Failed\n",i);
|
|
if(hamming(i,z)<=3)printf("Error: 0x%04X should map to 0x%04X\n",i,z);
|
|
}
|
|
else{
|
|
printf("0x%04X->0x%04X\n",i,y);
|
|
if(z!=y)printf("Error: 0x%04X should map to 0x%04X\n",i,z);
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
#endif
|