Algorithm-FEC
view release on metacpan or search on metacpan
/*
* i use malloc so many times, it is easier to put checks all in
* one place.
*/
static void *
my_malloc(int sz, char *err_string)
{
void *p = malloc( sz );
if (p == NULL) {
fprintf(stderr, "-- malloc failure allocating %s\n", err_string);
exit(1) ;
}
return p ;
}
#define NEW_GF_MATRIX(rows, cols) \
(gf *)my_malloc(rows * cols * sizeof(gf), " ## __LINE__ ## " )
/*
* initialize the data structures used for computations in GF.
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 ;
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);
return 0 ;
}
static int fec_initialized = 0 ;
static void init_fec()
{
TICK(ticks[0]);
generate_gf();
TOCK(ticks[0]);
DDB(fprintf(stderr, "generate_gf took %ldus\n", ticks[0]);)
TICK(ticks[0]);
init_mul_table();
TOCK(ticks[0]);
DDB(fprintf(stderr, "init_mul_table took %ldus\n", ticks[0]);)
fec_initialized = 1 ;
}
/*
* This section contains the proper FEC encoding/decoding routines.
* The encoding matrix is computed starting with a Vandermonde matrix,
* and then transforming it into a systematic matrix.
*/
struct fec_parms {
int k, n ; /* parameters of the code */
gf *enc_matrix ;
} ;
void
fec_free(struct fec_parms *p)
{
if (p==NULL) {
fprintf(stderr, "bad parameters to fec_free\n");
return ;
}
free(p->enc_matrix);
free(p);
}
/*
* create a new encoder, returning a descriptor. This contains k,n and
* the encoding matrix.
*/
{
int row, col ;
gf *p, *tmp_m ;
struct fec_parms *retval ;
if (fec_initialized == 0)
init_fec();
if (k > GF_SIZE + 1 || n > GF_SIZE + 1 || k > n ) {
fprintf(stderr, "Invalid parameters k %d n %d GF_SIZE %d\n",
k, n, GF_SIZE );
return NULL ;
}
retval = my_malloc(sizeof(struct fec_parms), "new_code");
retval->k = k ;
retval->n = n ;
retval->enc_matrix = NEW_GF_MATRIX(n, k);
tmp_m = NEW_GF_MATRIX(n, k);
/*
* fill the matrix with powers of field elements, starting from 0.
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".
*/
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
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 ;
}
#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 1.243 second using v1.01-cache-2.11-cpan-49f99fa48dc )