[Author Prev][Author Next][Thread Prev][Thread Next][Author Index][Thread Index]
[minion-cvs] Adding a copy of Luigi Rizzo"s software FEC code to con...
Update of /home/minion/cvsroot/src/minion/contrib/fec
In directory moria.mit.edu:/tmp/cvs-serv30511/fec
Added Files:
Makefile README fec.3 fec.S.980624a fec.S16.980624a fec.c
fec.h fec.s.980621e test.c
Log Message:
Adding a copy of Luigi Rizzo's software FEC code to contrib for K-of-N
fragmentation.
This code may eventually need to be replaced, because it performs poorly
with large N and K.
--- NEW FILE: Makefile ---
#
# makefile for vdm.
#
# fec.S.980624a is an optimized version for use on 486 and old pentium
# machines. It is only for GF_BITS=8 and generally does not work
# fast on systems with multiple instruction pipelines (PentiumPro,
# Pentium2)
# same for fec.S16.980624a , for use with GF_BITS=16
#
# gcc does something strange, so check the various opt. levels for
# best performance (or write addmul1 in assembly code).
#
# Standard compilation with -O9 works well for PentiumPro and Pentium2
# machines.
#
CC=gcc
# COPT= -O9 -funroll-loops -DGF_BITS=8
COPT= -O1 -DGF_BITS=8
CFLAGS=$(COPT) -Wall # -DTEST
SRCS= fec.c Makefile test.c fec.s.980621e \
fec.S.980624a \
fec.S16.980624a
DOCS= README fec.3
ALLSRCS= $(SRCS) $(DOCS) fec.h
fec: fec.o test.o
$(CC) $(CFLAGS) -o fec fec.o test.o
fec.o: fec.h fec.S
$(CC) $(CFLAGS) -c -o fec.o fec.S
fec.S: fec.c Makefile
$(CC) $(CFLAGS) -S -o fec.S fec.c
clean:
- rm *.core *.o fec.s fec
tgz: $(ALLSRCS)
tar cvzf vdm`date +%y%m%d`.tgz $(ALLSRCS)
--- NEW FILE: README ---
-------- Erasure codes based on Vandermonde matrices ---------
-------- (C) 1996-1998 Luigi Rizzo (luigi@iet.unipi.it) ---------
See fec.c for other contributors.
fec.c contains an implementation of an encoder/decoder for an erasure
code based on Vandermonde matrices computed over GF(2^m), m=2..16
PRINCIPLE OF OPERATION
The encoded data is computer as
y = E x
where x is a k-vector with source data, y is an n-vector with the
redundant info, and E is an n*k matrix derived from a Vandermonde
matrix. The code is systematic.
At the receiver, any subset y' of k elements from y allows the
reconstruction of the whole x by solving the system
y' = E' x
where E' is made of rows from E corresponding to the received
elements.
The complexity of matrix inversion is O(k*l^2) where l is the number
of elements not in x available at the receiver. This might seem
large, but data elements are in fact be packets of large size, so
the inversion cost can be amortized over the size of the packet.
For practical applications (k and l as large as 30, packet sizes
of 1KB) the cost can be neglected.
In addition, each of the l lost data packets has a reconstruction
cost O(k), (obviously) similar to the cost of the encoding phase.
Being the code systematic, you can express encoding and decoding costs
roughly as
Encoding speed: (c_e/l) MB/s
Decoding speed: (c_d/l) MB/s
PERFORMANCE
The values of c_d and c_e (similar) are very sensitive to issues
such as cache hit rate, memory speed, field size (8 or 16 bits),
etc. Also some machines are better than others in working with
bytes or 16-bit words. With the June'98 implementation I have
measured the following values for c_e and c_d (8-bit version; 16-bit
version has a penalty between 3 and 4).
Hardware C version Optimized version(*)
PentiumII 266 62 33
PentiumPRO 200 56 30
Pentium133 14.5 18
486dx2/66 4.05 4.55
(*) The 'optimized' version has some manual optimizations of the assembly
code generated by the compiler. It is generally faster for machines
with a single instruction pipeline, and generally slower for
machines with multiple pipelines.
See the manpage for detailed usage information.
--- NEW FILE: fec.3 ---
.\"
.\" Copyright 1998 by Luigi Rizzo, Dip. Ingegneria dell'Informazione,
.\" Universitaet Berlin. See the source code for copyright details.
.\" THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
.\"
.Dd July 15, 1998
.Dt FEC 3
.Os
.Sh NAME
.Nm fec_new, fec_encode, fec_encode, fec_free
.Nd An erasure code in GF(2^m)
.Sh SYNOPSIS
.Fd #include <fec.h>
.Ft void *
.Fn fec_new "int k" "int n"
.Ft void
.Fn fec_encode "void *code" "void *data[]" "void *dst" "int i" "int sz"
.Ft int
.Fn fec_decode "void *code" "void *data[]" "int i[]" "int sz"
.Ft void *
.Fn fec_free "void *code"
.Sh "DESCRIPTION"
This library implements a simple (n,k)
erasure code based on Vandermonde matrices.
The encoder takes
.Fa k
packets of size
.Fa sz
each, and is able to produce up to
.Fa n
different encoded packets, numbered from 0 to n-1,
such that any subset of
.Fa k
of them permits reconstruction of the original data.
.Pp
The data structures necessary for the encoding/decoding must
first be created using calling
.Fn fec_new
with the desired parameters. The code descriptor returned by the function
must be passed to other functions, and destroyed calling
.Fn fec_free
.Pp
Allowed values for k and n depend on a compile-time value
of
.Fa GF_BITS
and must be k <= n <= 2^GF_BITS.
Best performance is achieved with GF_BITS=8, although the code supports
also GF_BITS=16.
.Pp
Encoding is done by calling
.Fn fec_encode
and passing it pointers to the code descriptor, the source and
destination data packets, the index of the packet to be produced,
and the size of the packet.
.Pp Decoding is done calling
.Fn fec_decode
with pointers to the code, received packets, indexes of received
packets, and packet size. Decoding is done in place, possibly
shuffling the arrays passed as parameters. Decoding is deterministic
as long as the received packets are different. The decoding procedure
does some limited testing on this and returns if parameters are
invalid.
.Sh EXAMPLE
.nf
#include <fec.h>
/*
* example of sender code
*/
void *code ;
int n, k ;
void *src[] ;
void *pkt ;
code = new_code (k, n );
for (i = 0 ; i < k ; i++ )
src[i] = .. pointer to i-th source packet ..
for (each packet to transmit) {
i = ... index of the packet ;
fec_encode(code, src, pkt, i, size) ;
.. use packet in pkt
}
fec_free(code) ;
/*
* example of receiver code
*/
void *code ;
int n, k ;
void *data[] ;
int *ix[] ;
code = new_code (k, n );
for (i = 0 ; i < k ; i++ ) {
... receive a new packet ...
data[i] = .. pointer to i-th source packet ..
ix[i] = .. index of i-th source packet ..
}
fec_decode(code, data, ix, size) ;
/*
* now data[] has pointers to the source packets
*/
.SH BUGS
Please direct bug reports to luigi@iet.unipi.it .
.Sh "SEE ALSO"
--- NEW FILE: fec.S.980624a ---
(This appears to be a binary file; contents omitted.)
--- NEW FILE: fec.S16.980624a ---
(This appears to be a binary file; contents omitted.)
--- NEW FILE: fec.c ---
/*
* 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
*
* 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.
*/
/*
* The following parameter defines how many bits are used for
* field elements. The code supports any value from 2 to 16
* but fastest operation is achieved with 8 bit elements
* This is the only parameter you may want to change.
*/
#ifndef GF_BITS
#define GF_BITS 8 /* code over GF(2**GF_BITS) - change to suit */
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/*
* 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 unsigned char gf;
#else
typedef unsigned short 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;
}
#define SWAP(a,b,t) {t tmp; tmp=a; a=b; b=tmp;}
/*
* gf_mul(x,y) multiplies two numbers. If GF_BITS<=8, it is much
* faster to use a multiplication table.
*
* USE_GF_MULC, GF_MULC0(c) and GF_ADDMULC(x) can be used when multiplying
* many numbers by the same constant. In this case the first
* call sets the constant, and others perform the multiplications.
* A value related to the multiplication is held in a local variable
* declared with USE_GF_MULC . See usage in addmul1().
*/
#if (GF_BITS <= 8)
static gf gf_mul_table[GF_SIZE + 1][GF_SIZE + 1];
#define gf_mul(x,y) gf_mul_table[x][y]
#define USE_GF_MULC register gf * __gf_mulc_
#define GF_MULC0(c) __gf_mulc_ = gf_mul_table[c]
#define GF_ADDMULC(dst, x) dst ^= __gf_mulc_[x]
static void
init_mul_table()
{
int i, j;
for (i=0; i< GF_SIZE+1; i++)
for (j=0; j< GF_SIZE+1; j++)
gf_mul_table[i][j] = gf_exp[modnn(gf_log[i] + gf_log[j]) ] ;
for (j=0; j< GF_SIZE+1; j++)
gf_mul_table[0][j] = gf_mul_table[j][0] = 0;
}
#else /* GF_BITS > 8 */
static inline gf
gf_mul(x,y)
{
if ( (x) == 0 || (y)==0 ) return 0;
return gf_exp[gf_log[x] + gf_log[y] ] ;
}
#define init_mul_table()
#define USE_GF_MULC register gf * __gf_mulc_
#define GF_MULC0(c) __gf_mulc_ = &gf_exp[ gf_log[c] ]
#define GF_ADDMULC(dst, x) { if (x) dst ^= __gf_mulc_[ gf_log[x] ] ; }
#endif
/*
* Generate GF(2**m) from the irreducible polynomial p(X) in p[0]..p[m]
* Lookup tables:
* index->polynomial form gf_exp[] contains j= \alpha^i;
* polynomial form -> index form gf_log[ j = \alpha^i ] = i
* \alpha=x is the primitive element of GF(2^m)
*
* For efficiency, gf_exp[] has size 2*GF_SIZE, so that a simple
* multiplication of two numbers can be resolved without calling modnn
*/
/*
* 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.
*/
static void
generate_gf(void)
{
int i;
gf mask;
char *Pp = allPp[GF_BITS] ;
mask = 1; /* x ** 0 = 1 */
gf_exp[GF_BITS] = 0; /* will be updated at the end of the 1st loop */
/*
* first, generate the (polynomial representation of) powers of \alpha,
* which are stored in gf_exp[i] = \alpha ** i .
* At the same time build gf_log[gf_exp[i]] = i .
* The first GF_BITS powers are simply bits shifted to the left.
*/
for (i = 0; i < GF_BITS; i++, mask <<= 1 ) {
gf_exp[i] = mask;
gf_log[gf_exp[i]] = i;
/*
* If Pp[i] == 1 then \alpha ** i occurs in poly-repr
* gf_exp[GF_BITS] = \alpha ** GF_BITS
*/
if ( Pp[i] == '1' )
gf_exp[GF_BITS] ^= mask;
}
/*
* now gf_exp[GF_BITS] = \alpha ** GF_BITS is complete, so can als
* compute its inverse.
*/
gf_log[gf_exp[GF_BITS]] = GF_BITS;
/*
* Poly-repr of \alpha ** (i+1) is given by poly-repr of
* \alpha ** i shifted left one-bit and accounting for any
* \alpha ** GF_BITS term that may occur when poly-repr of
* \alpha ** i is shifted.
*/
mask = 1 << (GF_BITS - 1 ) ;
for (i = GF_BITS + 1; i < GF_SIZE; i++) {
if (gf_exp[i - 1] >= mask)
gf_exp[i] = gf_exp[GF_BITS] ^ ((gf_exp[i - 1] ^ mask) << 1);
else
gf_exp[i] = gf_exp[i - 1] << 1;
gf_log[gf_exp[i]] = i;
}
/*
* log(0) is not defined, so use a special value
*/
gf_log[0] = GF_SIZE ;
/* set the extended gf_exp values for fast multiply */
for (i = 0 ; i < GF_SIZE ; i++)
gf_exp[i + GF_SIZE] = gf_exp[i] ;
/*
* again special cases. 0 has no inverse. This used to
* be initialized to GF_SIZE, but it should make no difference
* since noone is supposed to read from here.
*/
inverse[0] = 0 ;
inverse[1] = 1;
for (i=2; i<=GF_SIZE; i++)
inverse[i] = gf_exp[GF_SIZE-gf_log[i]];
}
/*
* Various linear algebra operations that i use often.
*/
/*
* addmul() computes dst[] = dst[] + c * src[]
* This is used often, so better optimize it! Currently the loop is
* unrolled 16 times, a good value for 486 and pentium-class machines.
* The case c=0 is also optimized, whereas c=1 is not. These
* calls are unfrequent in my typical apps so I did not bother.
*
* Note that gcc on
*/
#define addmul(dst, src, c, sz) \
if (c != 0) addmul1(dst, src, c, sz)
#define UNROLL 16 /* 1, 4, 8, 16 */
static void
addmul1(gf *dst1, gf *src1, gf c, int sz)
{
USE_GF_MULC ;
register gf *dst = dst1, *src = src1 ;
gf *lim = &dst[sz - UNROLL + 1] ;
GF_MULC0(c) ;
#if (UNROLL > 1) /* unrolling by 8/16 is quite effective on the pentium */
for (; dst < lim ; dst += UNROLL, src += UNROLL ) {
GF_ADDMULC( dst[0] , src[0] );
GF_ADDMULC( dst[1] , src[1] );
GF_ADDMULC( dst[2] , src[2] );
GF_ADDMULC( dst[3] , src[3] );
#if (UNROLL > 4)
GF_ADDMULC( dst[4] , src[4] );
GF_ADDMULC( dst[5] , src[5] );
GF_ADDMULC( dst[6] , src[6] );
GF_ADDMULC( dst[7] , src[7] );
#endif
#if (UNROLL > 8)
GF_ADDMULC( dst[8] , src[8] );
GF_ADDMULC( dst[9] , src[9] );
GF_ADDMULC( dst[10] , src[10] );
GF_ADDMULC( dst[11] , src[11] );
GF_ADDMULC( dst[12] , src[12] );
GF_ADDMULC( dst[13] , src[13] );
GF_ADDMULC( dst[14] , src[14] );
GF_ADDMULC( dst[15] , src[15] );
#endif
}
#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)
*/
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] ;
}
for (col = 0 ; col < k ; col++ )
src[col*k + row] = gf_mul(inverse[t], b[col] );
}
free(c) ;
free(b) ;
free(p) ;
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.
*/
#define FEC_MAGIC 0xFECC0DEC
struct fec_parms {
u_long magic ;
int k, n ; /* parameters of the code */
gf *enc_matrix ;
} ;
void
fec_free(struct fec_parms *p)
{
if (p==NULL ||
p->magic != ( ( (FEC_MAGIC ^ p->k) ^ p->n) ^ (int)(p->enc_matrix)) ) {
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.
*/
struct fec_parms *
fec_new(int k, int n)
{
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);
retval->magic = ( ( FEC_MAGIC ^ k) ^ n) ^ (int)(retval->enc_matrix) ;
tmp_m = NEW_GF_MATRIX(n, k);
/*
* fill the matrix with powers of field elements, starting from 0.
* The first row is special, cannot be computed with exp. table.
*/
tmp_m[0] = 1 ;
for (col = 1; col < k ; col++)
tmp_m[col] = 0 ;
for (p = tmp_m + k, row = 0; row < n-1 ; row++, p += k) {
for ( col = 0 ; col < k ; col ++ )
p[col] = gf_exp[modnn(row*col)];
}
/*
* quick code to build systematic matrix: invert the top
* k*k vandermonde matrix, multiply right the bottom n-k rows
* by the inverse, and construct the identity matrix at the top.
*/
TICK(ticks[3]);
invert_vdm(tmp_m, k); /* much faster than invert_mat */
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]);
}
}
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 */
--- NEW FILE: fec.h ---
/*
* fec.c -- forward error correction based on Vandermonde matrices
* 980614
* (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
*
* 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.
*/
/*
* The following parameter defines how many bits are used for
* field elements. The code supports any value from 2 to 16
* but fastest operation is achieved with 8 bit elements
* This is the only parameter you may want to change.
*/
#ifndef GF_BITS
#define GF_BITS 8 /* code over GF(2**GF_BITS) - change to suit */
#endif
#define GF_SIZE ((1 << GF_BITS) - 1) /* powers of \alpha */
void fec_free(void *p) ;
void * fec_new(int k, int n) ;
void init_fec() ;
void fec_encode(void *code, void *src[], void *dst, int index, int sz) ;
int fec_decode(void *code, void *pkt[], int index[], int sz) ;
/* end of file */
--- NEW FILE: fec.s.980621e ---
(This appears to be a binary file; contents omitted.)
--- NEW FILE: test.c ---
/*
* test.c -- test code for FEC library
*
* (C) 1997-98 Luigi Rizzo (luigi@iet.unipi.it)
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "fec.h"
/*
* 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
*/
#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 */
void *
my_malloc(int sz, char *s)
{
void *p = malloc(sz) ;
if (p != NULL)
return p ;
fprintf(stderr, "test: malloc failure for %d bytes in <%s>\n",
sz, s);
exit(1);
}
int
pr_matrix(void *m1, int rows, int cols, char *s)
{
int r, c;
#if GF_BITS >8
u_char *m = m1 ;
#else
u_short *m = m1 ;
#endif
fprintf(stderr,"%s\n", s);
for (r=0; r<rows; r++) {
for (c=0; c<cols; c++)
fprintf(stderr,"%3d ",m[r*cols+c]);
fprintf(stderr,"\n");
}
fprintf(stderr,"\n");
return 0;
}
/*
* the following is only test code and can be safely omitted
* in applications.
* Creates k packets of size sz of random data, encodes them,
* and tries to decode.
* Index contains the permutation entry.
*/
int
test_decode(void *code, int k, int index[], int sz, char *s)
{
int errors;
int reconstruct = 0 ;
int item, i ;
static int prev_k = 0, prev_sz = 0;
static u_char **d_original = NULL, **d_src = NULL ;
if (sz < 1 || sz > 8192) {
fprintf(stderr, "test_decode: size %d invalid, must be 1..8K\n",
sz);
return 1 ;
}
if (k < 1 || k > GF_SIZE + 1) {
fprintf(stderr, "test_decode: k %d invalid, must be 1..%d\n",
k, GF_SIZE + 1 );
return 2 ;
}
if (prev_k != k || prev_sz != sz) {
if (d_original != NULL) {
for (i = 0 ; i < prev_k ; i++ ) {
free(d_original[i]);
free(d_src[i]);
}
free(d_original);
free(d_src);
d_original = NULL ;
d_src = NULL ;
}
}
prev_k = k ;
prev_sz = sz ;
if (d_original == NULL) {
d_original = my_malloc(k * sizeof(void *), "d_original ptr");
d_src = my_malloc(k * sizeof(void *), "d_src ptr");
for (i = 0 ; i < k ; i++ ) {
d_original[i] = my_malloc(sz, "d_original data");
d_src[i] = my_malloc(sz, "d_src data");
}
/*
* build sample data
*/
for (i = 0 ; i < k ; i++ ) {
for (item=0; item < sz; item++)
d_original[i][item] = ((item ^ i) + 3) & GF_SIZE;
}
}
errors = 0 ;
for( i = 0 ; i < k ; i++ )
if (index[i] >= k ) reconstruct ++ ;
TICK(ticks[2]);
for( i = 0 ; i < k ; i++ )
fec_encode(code, d_original, d_src[i], index[i], sz );
TOCK(ticks[2]);
TICK(ticks[1]);
if (fec_decode(code, d_src, index, sz)) {
fprintf(stderr, "detected singular matrix for %s \n", s);
return 1 ;
}
TOCK(ticks[1]);
for (i=0; i<k; i++)
if (bcmp(d_original[i], d_src[i], sz )) {
errors++;
fprintf(stderr, "error reconstructing block %d\n", i);
}
if (errors)
fprintf(stderr, "Errors reconstructing %d blocks out of %d\n",
errors, k);
fprintf(stderr,
" k %3d, l %3d c_enc %10.6f MB/s c_dec %10.6f MB/s \r",
k, reconstruct,
(double)(k * sz * reconstruct)/(double)ticks[2],
(double)(k * sz * reconstruct)/(double)ticks[1]);
return errors ;
}
#if 0
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
#define KK 64 /* 255 */
#define SZ 1024
int
main(int argc, char *argv[])
{
char buf[256];
void *code ;
int kk ;
int i ;
int *ixs ;
int lim = GF_SIZE + 1 ;
if (lim > 1024) lim = 1024 ;
#if 0
test_gf();
#endif
for ( kk = KK ; kk > 2 ; kk-- ) {
code = fec_new(kk, lim);
ixs = my_malloc(kk * sizeof(int), "ixs" );
for (i=0; i<kk; i++) ixs[i] = kk - i ;
sprintf(buf, "kk=%d, kk - i", kk);
test_decode(code, kk, ixs, SZ, buf);
for (i=0; i<kk; i++) ixs[i] = i ;
test_decode(code, kk, ixs, SZ, "i");
if (0) {
for (i=0; i<kk; i++) ixs[i] = i ;
ixs[0] = ixs[kk/2] ;
test_decode(code, kk, ixs, SZ, "0 = 1 (error expected)");
}
if (0)
for (i= lim-1 ; i >= kk ; i--) {
int j ;
for (j=0; j<KK; j++) ixs[j] = kk - j ;
ixs[0] = i ;
test_decode(code, kk, ixs, SZ, "0 = big");
}
if (0)
for (i= lim - kk ; i >= 0 && i>= lim - kk - 4 ; i--) {
int j ;
for (j=0; j<kk; j++)
ixs[j] = kk -1 - j + i ;
test_decode(code, kk, ixs, SZ, "shifted j");
}
if (1) {
int j, max_i0 = KK/2 ;
if (max_i0 + KK > lim)
max_i0 = lim - KK ;
for (i= 0 ; i <= max_i0 ; i++) {
for (j=0; j<kk; j++)
ixs[j] = j + i ;
test_decode(code, kk, ixs, SZ, "shifted j");
}
}
fprintf(stderr, "\n");
free(ixs);
fec_free(code);
}
return 0;
}