Algorithm-FEC
view release on metacpan or search on metacpan
/*
* fec.c -- forward error correction based on Vandermonde matrices
* 980624
* (C) 1997-98 Luigi Rizzo (luigi@iet.unipi.it)
*
* Portions derived from code by Phil Karn (karn@ka9q.ampr.org),
* Robert Morelos-Zaragoza (robert@spectra.eng.hawaii.edu) and Hari
* Thirumoorthy (harit@spectra.eng.hawaii.edu), Aug 1995
* modified by Marc Lehmann <fec@schmorp.de>, Sep 2003.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
*
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above
* copyright notice, this list of conditions and the following
* disclaimer in the documentation and/or other materials
* provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
* THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
* PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS
* BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
* OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
* OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
* TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
* OF SUCH DAMAGE.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#define MSDOS /* LEAVE THIS IN PLACE EVEN ON UNIX! */
/*
* compatibility stuff
*/
#ifdef MSDOS /* but also for others, e.g. sun... */
#define NEED_BCOPY
#define bcmp(a,b,n) memcmp(a,b,n)
#endif
#ifdef NEED_BCOPY
#define bcopy(s, d, siz) memcpy((d), (s), (siz))
#define bzero(d, siz) memset((d), '\0', (siz))
#endif
/*
* stuff used for testing purposes only
*/
#ifdef TEST
#define DEB(x)
#define DDB(x) x
#define DEBUG 0 /* minimal debugging */
#ifdef MSDOS
#include <time.h>
struct timeval {
unsigned long ticks;
};
#define gettimeofday(x, dummy) { (x)->ticks = clock() ; }
#define DIFF_T(a,b) (1+ 1000000*(a.ticks - b.ticks) / CLOCKS_PER_SEC )
typedef unsigned long u_long ;
typedef unsigned short u_short ;
#else /* typically, unix systems */
#include <sys/time.h>
#define DIFF_T(a,b) \
(1+ 1000000*(a.tv_sec - b.tv_sec) + (a.tv_usec - b.tv_usec) )
#endif
#define TICK(t) \
{struct timeval x ; \
gettimeofday(&x, NULL) ; \
t = x.tv_usec + 1000000* (x.tv_sec & 0xff ) ; \
}
#define TOCK(t) \
{ u_long t1 ; TICK(t1) ; \
if (t1 < t) t = 256000000 + t1 - t ; \
else t = t1 - t ; \
if (t == 0) t = 1 ;}
u_long ticks[10]; /* vars for timekeeping */
#else
#define DEB(x)
#define DDB(x)
#define TICK(x)
#define TOCK(x)
#endif /* TEST */
/*
* You should not need to change anything beyond this point.
* The first part of the file implements linear algebra in GF.
*
* gf is the type used to store an element of the Galois Field.
* Must constain at least GF_BITS bits.
*
* Note: unsigned char will work up to GF(256) but int seems to run
* faster on the Pentium. We use int whenever have to deal with an
* index, since they are generally faster.
*/
#if (GF_BITS < 2 && GF_BITS >16)
#error "GF_BITS must be 2 .. 16"
#endif
#if (GF_BITS <= 8)
typedef uint8_t gf;
#else
typedef uint16_t gf;
#endif
#define GF_SIZE ((1 << GF_BITS) - 1) /* powers of \alpha */
/*
* Primitive polynomials - see Lin & Costello, Appendix A,
* and Lee & Messerschmitt, p. 453.
*/
static char *allPp[] = { /* GF_BITS polynomial */
NULL, /* 0 no code */
NULL, /* 1 no code */
"111", /* 2 1+x+x^2 */
"1101", /* 3 1+x+x^3 */
"11001", /* 4 1+x+x^4 */
"101001", /* 5 1+x^2+x^5 */
"1100001", /* 6 1+x+x^6 */
"10010001", /* 7 1 + x^3 + x^7 */
"101110001", /* 8 1+x^2+x^3+x^4+x^8 */
"1000100001", /* 9 1+x^4+x^9 */
"10010000001", /* 10 1+x^3+x^10 */
"101000000001", /* 11 1+x^2+x^11 */
"1100101000001", /* 12 1+x+x^4+x^6+x^12 */
"11011000000001", /* 13 1+x+x^3+x^4+x^13 */
"110000100010001", /* 14 1+x+x^6+x^10+x^14 */
"1100000000000001", /* 15 1+x+x^15 */
"11010000000010001" /* 16 1+x+x^3+x^12+x^16 */
};
/*
* To speed up computations, we have tables for logarithm, exponent
* and inverse of a number. If GF_BITS <= 8, we use a table for
* multiplication as well (it takes 64K, no big deal even on a PDA,
* especially because it can be pre-initialized an put into a ROM!),
* otherwhise we use a table of logarithms.
* In any case the macro gf_mul(x,y) takes care of multiplications.
*/
static gf gf_exp[2*GF_SIZE]; /* index->poly form conversion table */
static int gf_log[GF_SIZE + 1]; /* Poly->index form conversion table */
static gf inverse[GF_SIZE+1]; /* inverse of field elem. */
/* inv[\alpha**i]=\alpha**(GF_SIZE-i-1) */
/*
* modnn(x) computes x % GF_SIZE, where GF_SIZE is 2**GF_BITS - 1,
* without a slow divide.
*/
static inline gf
modnn(int x)
{
while (x >= GF_SIZE) {
x -= GF_SIZE;
x = (x >> GF_BITS) + (x & GF_SIZE);
}
return x;
}
#endif
lim += UNROLL - 1 ;
for (; dst < lim; dst++, src++ ) /* final components */
GF_ADDMULC( *dst , *src );
}
/*
* computes C = AB where A is n*k, B is k*m, C is n*m
*/
static void
matmul(gf *a, gf *b, gf *c, int n, int k, int m)
{
int row, col, i ;
for (row = 0; row < n ; row++) {
for (col = 0; col < m ; col++) {
gf *pa = &a[ row * k ];
gf *pb = &b[ col ];
gf acc = 0 ;
for (i = 0; i < k ; i++, pa++, pb += m )
acc ^= gf_mul( *pa, *pb ) ;
c[ row * m + col ] = acc ;
}
}
}
#ifdef DEBUG
/*
* returns 1 if the square matrix is identiy
* (only for test)
*/
static int
is_identity(gf *m, int k)
{
int row, col ;
for (row=0; row<k; row++)
for (col=0; col<k; col++)
if ( (row==col && *m != 1) ||
(row!=col && *m != 0) )
return 0 ;
else
m++ ;
return 1 ;
}
#endif /* debug */
/*
* invert_mat() takes a matrix and produces its inverse
* k is the size of the matrix.
* (Gauss-Jordan, adapted from Numerical Recipes in C)
* Return non-zero if singular.
*/
DEB( int pivloops=0; int pivswaps=0 ; /* diagnostic */)
static int
invert_mat(gf *src, int k)
{
gf c, *p ;
int irow, icol, row, col, i, ix ;
int error = 1 ;
int *indxc = my_malloc(k*sizeof(int), "indxc");
int *indxr = my_malloc(k*sizeof(int), "indxr");
int *ipiv = my_malloc(k*sizeof(int), "ipiv");
gf *id_row = NEW_GF_MATRIX(1, k);
gf *temp_row = NEW_GF_MATRIX(1, k);
bzero(id_row, k*sizeof(gf));
DEB( pivloops=0; pivswaps=0 ; /* diagnostic */ )
/*
* ipiv marks elements already used as pivots.
*/
for (i = 0; i < k ; i++)
ipiv[i] = 0 ;
for (col = 0; col < k ; col++) {
gf *pivot_row ;
/*
* Zeroing column 'col', look for a non-zero element.
* First try on the diagonal, if it fails, look elsewhere.
*/
irow = icol = -1 ;
if (ipiv[col] != 1 && src[col*k + col] != 0) {
irow = col ;
icol = col ;
goto found_piv ;
}
for (row = 0 ; row < k ; row++) {
if (ipiv[row] != 1) {
for (ix = 0 ; ix < k ; ix++) {
DEB( pivloops++ ; )
if (ipiv[ix] == 0) {
if (src[row*k + ix] != 0) {
irow = row ;
icol = ix ;
goto found_piv ;
}
} else if (ipiv[ix] > 1) {
fprintf(stderr, "singular matrix\n");
goto fail ;
}
}
}
}
if (icol == -1) {
fprintf(stderr, "XXX pivot not found!\n");
goto fail ;
}
found_piv:
++(ipiv[icol]) ;
/*
* swap rows irow and icol, so afterwards the diagonal
* element will be correct. Rarely done, not worth
* optimizing.
*/
if (irow != icol) {
for (ix = 0 ; ix < k ; ix++ ) {
SWAP( src[irow*k + ix], src[icol*k + ix], gf) ;
}
}
indxr[col] = irow ;
indxc[col] = icol ;
pivot_row = &src[icol*k] ;
c = pivot_row[icol] ;
if (c == 0) {
fprintf(stderr, "singular matrix 2\n");
goto fail ;
}
if (c != 1 ) { /* otherwhise this is a NOP */
/*
* this is done often , but optimizing is not so
* fruitful, at least in the obvious ways (unrolling)
*/
DEB( pivswaps++ ; )
c = inverse[ c ] ;
pivot_row[icol] = 1 ;
for (ix = 0 ; ix < k ; ix++ )
pivot_row[ix] = gf_mul(c, pivot_row[ix] );
}
/*
* from all rows, remove multiples of the selected row
* to zero the relevant entry (in fact, the entry is not zero
* because we know it must be zero).
* (Here, if we know that the pivot_row is the identity,
* we can optimize the addmul).
*/
id_row[icol] = 1;
if (bcmp(pivot_row, id_row, k*sizeof(gf)) != 0) {
for (p = src, ix = 0 ; ix < k ; ix++, p += k ) {
if (ix != icol) {
c = p[icol] ;
p[icol] = 0 ;
addmul(p, pivot_row, c, k );
}
}
}
id_row[icol] = 0;
} /* done all columns */
for (col = k-1 ; col >= 0 ; col-- ) {
if (indxr[col] <0 || indxr[col] >= k)
fprintf(stderr, "AARGH, indxr[col] %d\n", indxr[col]);
else if (indxc[col] <0 || indxc[col] >= k)
fprintf(stderr, "AARGH, indxc[col] %d\n", indxc[col]);
else
if (indxr[col] != indxc[col] ) {
for (row = 0 ; row < k ; row++ ) {
SWAP( src[row*k + indxr[col]], src[row*k + indxc[col]], gf) ;
}
}
}
error = 0 ;
fail:
free(indxc);
free(indxr);
free(ipiv);
free(id_row);
free(temp_row);
return error ;
}
/*
* fast code for inverting a vandermonde matrix.
* XXX NOTE: It assumes that the matrix
* is not singular and _IS_ a vandermonde matrix. Only uses
* the second column of the matrix, containing the p_i's.
*
* Algorithm borrowed from "Numerical recipes in C" -- sec.2.8, but
* largely revised for my purposes.
* p = coefficients of the matrix (p_i)
* q = values of the polynomial (known)
*/
static int
invert_vdm(gf *src, int k)
{
int i, j, row, col ;
gf *b, *c, *p;
gf t, xx ;
if (k == 1) /* degenerate case, matrix must be p^0 = 1 */
return 0 ;
/*
* c holds the coefficient of P(x) = Prod (x - p_i), i=0..k-1
* b holds the coefficient for the matrix inversion
*/
c = NEW_GF_MATRIX(1, k);
b = NEW_GF_MATRIX(1, k);
p = NEW_GF_MATRIX(1, k);
for ( j=1, i = 0 ; i < k ; i++, j+=k ) {
c[i] = 0 ;
p[i] = src[j] ; /* p[i] */
}
/*
* construct coeffs. recursively. We know c[k] = 1 (implicit)
* and start P_0 = x - p_0, then at each stage multiply by
* x - p_i generating P_i = x P_{i-1} - p_i P_{i-1}
* After k steps we are done.
*/
c[k-1] = p[0] ; /* really -p(0), but x = -x in GF(2^m) */
for (i = 1 ; i < k ; i++ ) {
gf p_i = p[i] ; /* see above comment */
for (j = k-1 - ( i - 1 ) ; j < k-1 ; j++ )
c[j] ^= gf_mul( p_i, c[j+1] ) ;
c[k-1] ^= p_i ;
}
for (row = 0 ; row < k ; row++ ) {
/*
* synthetic division etc.
*/
xx = p[row] ;
t = 1 ;
b[k-1] = 1 ; /* this is in fact c[k] */
for (i = k-2 ; i >= 0 ; i-- ) {
b[i] = c[i+1] ^ gf_mul(xx, b[i+1]) ;
t = gf_mul(xx, t) ^ b[i] ;
matmul(tmp_m + k*k, tmp_m, retval->enc_matrix + k*k, n - k, k, k);
/*
* the upper matrix is I so do not bother with a slow multiply
*/
bzero(retval->enc_matrix, k*k*sizeof(gf) );
for (p = retval->enc_matrix, col = 0 ; col < k ; col++, p += k+1 )
*p = 1 ;
free(tmp_m);
TOCK(ticks[3]);
DDB(fprintf(stderr, "--- %ld us to build encoding matrix\n",
ticks[3]);)
DEB(pr_matrix(retval->enc_matrix, n, k, "encoding_matrix");)
return retval ;
}
/*
* fec_encode accepts as input pointers to n data packets of size sz,
* and produces as output a packet pointed to by fec, computed
* with index "index".
*/
void
fec_encode(struct fec_parms *code, gf *src[], gf *fec, int index, int sz)
{
int i, k = code->k ;
gf *p ;
if (GF_BITS > 8)
sz /= 2 ;
if (index < k)
bcopy(src[index], fec, sz*sizeof(gf) ) ;
else if (index < code->n) {
p = &(code->enc_matrix[index*k] );
bzero(fec, sz*sizeof(gf));
for (i = 0; i < k ; i++)
addmul(fec, src[i], p[i], sz ) ;
} else
fprintf(stderr, "Invalid index %d (max %d)\n",
index, code->n - 1 );
}
/*
* shuffle move src packets in their position
*/
static int
shuffle(gf *pkt[], int index[], int k)
{
int i;
for ( i = 0 ; i < k ; ) {
if (index[i] >= k || index[i] == i)
i++ ;
else {
/*
* put pkt in the right position (first check for conflicts).
*/
int c = index[i] ;
if (index[c] == c) {
DEB(fprintf(stderr, "\nshuffle, error at %d\n", i);)
return 1 ;
}
SWAP(index[i], index[c], int) ;
SWAP(pkt[i], pkt[c], gf *) ;
}
}
DEB( /* just test that it works... */
for ( i = 0 ; i < k ; i++ ) {
if (index[i] < k && index[i] != i) {
fprintf(stderr, "shuffle: after\n");
for (i=0; i<k ; i++) fprintf(stderr, "%3d ", index[i]);
fprintf(stderr, "\n");
return 1 ;
}
}
)
return 0 ;
}
/*
* build_decode_matrix constructs the encoding matrix given the
* indexes. The matrix must be already allocated as
* a vector of k*k elements, in row-major order
*/
static gf *
build_decode_matrix(struct fec_parms *code, gf *pkt[], int index[])
{
int i , k = code->k ;
gf *p, *matrix = NEW_GF_MATRIX(k, k);
TICK(ticks[9]);
for (i = 0, p = matrix ; i < k ; i++, p += k ) {
#if 1 /* this is simply an optimization, not very useful indeed */
if (index[i] < k) {
bzero(p, k*sizeof(gf) );
p[i] = 1 ;
} else
#endif
if (index[i] < code->n )
bcopy( &(code->enc_matrix[index[i]*k]), p, k*sizeof(gf) );
else {
fprintf(stderr, "decode: invalid index %d (max %d)\n",
index[i], code->n - 1 );
free(matrix) ;
return NULL ;
}
}
TICK(ticks[9]);
if (invert_mat(matrix, k)) {
free(matrix);
matrix = NULL ;
}
TOCK(ticks[9]);
return matrix ;
}
/*
* fec_decode receives as input a vector of packets, the indexes of
* packets, and produces the correct vector as output.
*
* Input:
* code: pointer to code descriptor
* pkt: pointers to received packets. They are modified
* to store the output packets (in place)
* index: pointer to packet indexes (modified)
* sz: size of each packet
*/
int
fec_decode(struct fec_parms *code, gf *pkt[], int index[], int sz)
{
gf *m_dec ;
gf **new_pkt ;
int row, col , k = code->k ;
if (GF_BITS > 8)
sz /= 2 ;
if (shuffle(pkt, index, k)) /* error if true */
return 1 ;
m_dec = build_decode_matrix(code, pkt, index);
if (m_dec == NULL)
return 1 ; /* error */
/*
* do the actual decoding
*/
new_pkt = my_malloc (k * sizeof (gf * ), "new pkt pointers" );
for (row = 0 ; row < k ; row++ ) {
if (index[row] >= k) {
new_pkt[row] = my_malloc (sz * sizeof (gf), "new pkt buffer" );
bzero(new_pkt[row], sz * sizeof(gf) ) ;
for (col = 0 ; col < k ; col++ )
addmul(new_pkt[row], pkt[col], m_dec[row*k + col], sz) ;
}
}
/*
* move pkts to their final destination
*/
for (row = 0 ; row < k ; row++ ) {
if (index[row] >= k) {
bcopy(new_pkt[row], pkt[row], sz*sizeof(gf));
free(new_pkt[row]);
index[row] = row;
}
}
free(new_pkt);
free(m_dec);
return 0;
}
/*********** end of FEC code -- beginning of test code ************/
#if (TEST || DEBUG)
void
test_gf()
{
int i ;
/*
* test gf tables. Sufficiently tested...
*/
for (i=0; i<= GF_SIZE; i++) {
if (gf_exp[gf_log[i]] != i)
fprintf(stderr, "bad exp/log i %d log %d exp(log) %d\n",
i, gf_log[i], gf_exp[gf_log[i]]);
if (i != 0 && gf_mul(i, inverse[i]) != 1)
fprintf(stderr, "bad mul/inv i %d inv %d i*inv(i) %d\n",
i, inverse[i], gf_mul(i, inverse[i]) );
if (gf_mul(0,i) != 0)
fprintf(stderr, "bad mul table 0,%d\n",i);
if (gf_mul(i,0) != 0)
fprintf(stderr, "bad mul table %d,0\n",i);
}
}
#endif /* TEST */
( run in 0.457 second using v1.01-cache-2.11-cpan-f6376fbd888 )