fec.c

Go to the documentation of this file.
00001 
00035 #define FEC                     /* VR: added */
00036 
00037 #ifdef FEC
00038 
00039 #include "defines.h"
00040 #include "fec.h"
00041 
00042 #include <stdio.h>
00043 #include <stdlib.h>
00044 #include <string.h>
00045 
00046 #ifndef _MSC_VER
00047 #include <strings.h>    /* VR: added */
00048 #endif
00049 
00050 #include <sys/types.h>  /* VR: added */
00051 
00052 /*
00053  * compatibility stuff
00054  */
00055 #if (defined(MSDOS) || defined(_MSC_VER))       /* but also for others, e.g. sun... */
00056 #define NEED_BCOPY
00057 #define bcmp(a,b,n) memcmp(a,b,n)
00058 #endif
00059 
00060 #ifdef NEED_BCOPY
00061 #define bcopy(s, d, siz)        memcpy((d), (s), (siz))
00062 #define bzero(d, siz)   memset((d), '\0', (siz))
00063 #endif
00064 
00065 #ifndef u_long
00066         #define u_long unsigned long
00067 #endif
00068 
00069 /*
00070  * stuff used for testing purposes only
00071  */
00072 
00073 #ifdef TICK             /* VR: avoid a warning under Solaris */
00074 #undef TICK
00075 #endif
00076 
00077 #ifdef  TEST
00078 #define DEB(x)
00079 #define DDB(x) x
00080 #define DEBUG   4       /* minimal debugging */
00081 #if (defined(MSDOS) || defined(_MSC_VER))
00082 #include <time.h>
00083 struct timeval {
00084     unsigned long ticks;
00085 };
00086 #define gettimeofday(x, dummy) { (x)->ticks = clock() ; }
00087 #define DIFF_T(a,b) (1+ 1000000*(a.ticks - b.ticks) / CLOCKS_PER_SEC )
00088 typedef unsigned long u_long ;
00089 typedef unsigned short u_short ;
00090 #else /* typically, unix systems */
00091 #include <sys/time.h>
00092 #define DIFF_T(a,b) \
00093         (1+ 1000000*(a.tv_sec - b.tv_sec) + (a.tv_usec - b.tv_usec) )
00094 #endif
00095 
00096 #define TICK(t) \
00097         {struct timeval x ; \
00098         gettimeofday(&x, NULL) ; \
00099         t = x.tv_usec + 1000000* (x.tv_sec & 0xff ) ; \
00100         }
00101 #define TOCK(t) \
00102         { u_long t1 ; TICK(t1) ; \
00103           if (t1 < t) t = 256000000 + t1 - t ; \
00104           else t = t1 - t ; \
00105           if (t == 0) t = 1 ;}
00106 
00107 u_long ticks[10];       /* vars for timekeeping */
00108 #else
00109 #define DEB(x)
00110 #define DDB(x)
00111 #define TICK(x)
00112 #define TOCK(x)
00113 #endif /* TEST */
00114 
00115 /*
00116  * You should not need to change anything beyond this point.
00117  * The first part of the file implements linear algebra in GF.
00118  *
00119  * gf is the type used to store an element of the Galois Field.
00120  * Must constain at least GF_BITS bits.
00121  *
00122  * Note: unsigned char will work up to GF(256) but int seems to run
00123  * faster on the Pentium. We use int whenever have to deal with an
00124  * index, since they are generally faster.
00125  */
00126 #if (GF_BITS < 2  || GF_BITS >16)
00127 #error "GF_BITS must be 2 .. 16"
00128 #endif
00129 
00130 #if (GF_BITS <= 8)
00131 typedef unsigned char gf;
00132 #else
00133 typedef unsigned short gf;
00134 #endif
00135 
00136 // defined in fec.h
00137 //#define       GF_SIZE ((1 << GF_BITS) - 1)    /* powers of \alpha */
00138 
00139 /*
00140  * Primitive polynomials - see Lin & Costello, Appendix A,
00141  * and  Lee & Messerschmitt, p. 453.
00142  */
00143 static char *allPp[] = {    /* GF_BITS  polynomial              */
00144     NULL,                   /*  0       no code                 */
00145     NULL,                   /*  1       no code                 */
00146     "111",                  /*  2       1+x+x^2                 */
00147     "1101",                 /*  3       1+x+x^3                 */
00148     "11001",                /*  4       1+x+x^4                 */
00149     "101001",               /*  5       1+x^2+x^5               */
00150     "1100001",              /*  6       1+x+x^6                 */
00151     "10010001",             /*  7       1 + x^3 + x^7           */
00152     "101110001",            /*  8       1+x^2+x^3+x^4+x^8       */
00153     "1000100001",           /*  9       1+x^4+x^9               */
00154     "10010000001",          /* 10       1+x^3+x^10              */
00155     "101000000001",         /* 11       1+x^2+x^11              */
00156     "1100101000001",        /* 12       1+x+x^4+x^6+x^12        */
00157     "11011000000001",       /* 13       1+x+x^3+x^4+x^13        */
00158     "110000100010001",      /* 14       1+x+x^6+x^10+x^14       */
00159     "1100000000000001",     /* 15       1+x+x^15                */
00160     "11010000000010001"     /* 16       1+x+x^3+x^12+x^16       */
00161 };
00162 
00163 
00164 /*
00165  * To speed up computations, we have tables for logarithm, exponent
00166  * and inverse of a number. If GF_BITS <= 8, we use a table for
00167  * multiplication as well (it takes 64K, no big deal even on a PDA,
00168  * especially because it can be pre-initialized an put into a ROM!),
00169  * otherwhise we use a table of logarithms.
00170  * In any case the macro gf_mul(x,y) takes care of multiplications.
00171  */
00172 
00173 static gf gf_exp[2*GF_SIZE];    /* index->poly form conversion table    */
00174 static int gf_log[GF_SIZE + 1]; /* Poly->index form conversion table    */
00175 static gf inverse[GF_SIZE+1];   /* inverse of field elem.               */
00176                                 /* inv[\alpha**i]=\alpha**(GF_SIZE-i-1) */
00177 
00178 /*
00179  * modnn(x) computes x % GF_SIZE, where GF_SIZE is 2**GF_BITS - 1,
00180  * without a slow divide.
00181  */
00182 static gf
00183 modnn(int x)
00184 {
00185     while (x >= GF_SIZE) {
00186         x -= GF_SIZE;
00187         x = (x >> GF_BITS) + (x & GF_SIZE);
00188     }
00189     return x;
00190 }
00191 
00192 #define SWAP(a,b,t) {t tmp; tmp=a; a=b; b=tmp;}
00193 
00194 /*
00195  * gf_mul(x,y) multiplies two numbers. If GF_BITS<=8, it is much
00196  * faster to use a multiplication table.
00197  *
00198  * USE_GF_MULC, GF_MULC0(c) and GF_ADDMULC(x) can be used when multiplying
00199  * many numbers by the same constant. In this case the first
00200  * call sets the constant, and others perform the multiplications.
00201  * A value related to the multiplication is held in a local variable
00202  * declared with USE_GF_MULC . See usage in addmul1().
00203  */
00204 #if (GF_BITS <= 8)
00205 static gf gf_mul_table[GF_SIZE + 1][GF_SIZE + 1];
00206 
00207 #define gf_mul(x,y) gf_mul_table[x][y]
00208 
00209 #define USE_GF_MULC register gf * __gf_mulc_
00210 #define GF_MULC0(c) __gf_mulc_ = gf_mul_table[c]
00211 #define GF_ADDMULC(dst, x) dst ^= __gf_mulc_[x]
00212 
00213 static void
00214 init_mul_table()
00215 {
00216     int i, j;
00217     for (i=0; i< GF_SIZE+1; i++)
00218         for (j=0; j< GF_SIZE+1; j++)
00219             gf_mul_table[i][j] = gf_exp[modnn(gf_log[i] + gf_log[j]) ] ;
00220 
00221     for (j=0; j< GF_SIZE+1; j++)
00222             gf_mul_table[0][j] = gf_mul_table[j][0] = 0;
00223 }
00224 #else   /* GF_BITS > 8 */
00225 static inline gf
00226 gf_mul(x,y)
00227 {
00228     if ( (x) == 0 || (y)==0 ) return 0;
00229 
00230     return gf_exp[gf_log[x] + gf_log[y] ] ;
00231 }
00232 #define init_mul_table()
00233 
00234 #define USE_GF_MULC register gf * __gf_mulc_
00235 #define GF_MULC0(c) __gf_mulc_ = &gf_exp[ gf_log[c] ]
00236 #define GF_ADDMULC(dst, x) { if (x) dst ^= __gf_mulc_[ gf_log[x] ] ; }
00237 #endif
00238 
00239 /*
00240  * Generate GF(2**m) from the irreducible polynomial p(X) in p[0]..p[m]
00241  * Lookup tables:
00242  *     index->polynomial form           gf_exp[] contains j= \alpha^i;
00243  *     polynomial form -> index form    gf_log[ j = \alpha^i ] = i
00244  * \alpha=x is the primitive element of GF(2^m)
00245  *
00246  * For efficiency, gf_exp[] has size 2*GF_SIZE, so that a simple
00247  * multiplication of two numbers can be resolved without calling modnn
00248  */
00249 
00250 /*
00251  * i use malloc so many times, it is easier to put checks all in
00252  * one place.
00253  */
00254 static void *
00255 my_malloc(int sz, char *err_string)
00256 {
00257     void *p = malloc( sz );
00258     if (p == NULL) {
00259         fprintf( stderr, "-- malloc failure allocation %s\n", err_string );
00260         exit(1) ;
00261     }
00262     return p ;
00263 }
00264 
00265 #define NEW_GF_MATRIX(rows, cols) \
00266     (gf *)my_malloc(rows * cols * sizeof(gf), " ## __LINE__ ## " )
00267 
00268 /*
00269  * initialize the data structures used for computations in GF.
00270  */
00271 static void
00272 generate_gf(void)
00273 {
00274     int i;
00275     gf mask;
00276     char *Pp =  allPp[GF_BITS] ;
00277 
00278     mask = 1;   /* x ** 0 = 1 */
00279     gf_exp[GF_BITS] = 0; /* will be updated at the end of the 1st loop */
00280     /*
00281      * first, generate the (polynomial representation of) powers of \alpha,
00282      * which are stored in gf_exp[i] = \alpha ** i .
00283      * At the same time build gf_log[gf_exp[i]] = i .
00284      * The first GF_BITS powers are simply bits shifted to the left.
00285      */
00286     for (i = 0; i < GF_BITS; i++, mask <<= 1 ) {
00287         gf_exp[i] = mask;
00288         gf_log[gf_exp[i]] = i;
00289         /*
00290          * If Pp[i] == 1 then \alpha ** i occurs in poly-repr
00291          * gf_exp[GF_BITS] = \alpha ** GF_BITS
00292          */
00293         if ( Pp[i] == '1' )
00294             gf_exp[GF_BITS] ^= mask;
00295     }
00296     /*
00297      * now gf_exp[GF_BITS] = \alpha ** GF_BITS is complete, so can als
00298      * compute its inverse.
00299      */
00300     gf_log[gf_exp[GF_BITS]] = GF_BITS;
00301     /*
00302      * Poly-repr of \alpha ** (i+1) is given by poly-repr of
00303      * \alpha ** i shifted left one-bit and accounting for any
00304      * \alpha ** GF_BITS term that may occur when poly-repr of
00305      * \alpha ** i is shifted.
00306      */
00307     mask = 1 << (GF_BITS - 1 ) ;
00308     for (i = GF_BITS + 1; i < GF_SIZE; i++) {
00309         if (gf_exp[i - 1] >= mask)
00310             gf_exp[i] = gf_exp[GF_BITS] ^ ((gf_exp[i - 1] ^ mask) << 1);
00311         else
00312             gf_exp[i] = gf_exp[i - 1] << 1;
00313         gf_log[gf_exp[i]] = i;
00314     }
00315     /*
00316      * log(0) is not defined, so use a special value
00317      */
00318     gf_log[0] = GF_SIZE ;
00319     /* set the extended gf_exp values for fast multiply */
00320     for (i = 0 ; i < GF_SIZE ; i++)
00321         gf_exp[i + GF_SIZE] = gf_exp[i] ;
00322 
00323     /*
00324      * again special cases. 0 has no inverse. This used to
00325      * be initialized to GF_SIZE, but it should make no difference
00326      * since noone is supposed to read from here.
00327     */
00328     inverse[0] = 0 ;
00329     inverse[1] = 1;
00330     for (i=2; i<=GF_SIZE; i++)
00331         inverse[i] = gf_exp[GF_SIZE-gf_log[i]];
00332 }
00333 
00334 /*
00335  * Various linear algebra operations that i use often.
00336  */
00337 
00338 /*
00339  * addmul() computes dst[] = dst[] + c * src[]
00340  * This is used often, so better optimize it! Currently the loop is
00341  * unrolled 16 times, a good value for 486 and pentium-class machines.
00342  * The case c=0 is also optimized, whereas c=1 is not. These
00343  * calls are unfrequent in my typical apps so I did not bother.
00344  *
00345  * Note that gcc on
00346  */
00347 #define addmul(dst, src, c, sz) \
00348     if (c != 0) addmul1(dst, src, c, sz)
00349 
00350 #define UNROLL 16 /* 1, 4, 8, 16 */
00351 static void
00352 addmul1(gf *dst1, gf *src1, gf c, int sz)
00353 {
00354     USE_GF_MULC ;
00355     register gf *dst = dst1, *src = src1 ;
00356     gf *lim = &dst[sz - UNROLL + 1] ;
00357 
00358     GF_MULC0(c) ;
00359 
00360 #if (UNROLL > 1) /* unrolling by 8/16 is quite effective on the pentium */
00361     for (; dst < lim ; dst += UNROLL, src += UNROLL ) {
00362         GF_ADDMULC( dst[0] , src[0] );
00363         GF_ADDMULC( dst[1] , src[1] );
00364         GF_ADDMULC( dst[2] , src[2] );
00365         GF_ADDMULC( dst[3] , src[3] );
00366 #if (UNROLL > 4)
00367         GF_ADDMULC( dst[4] , src[4] );
00368         GF_ADDMULC( dst[5] , src[5] );
00369         GF_ADDMULC( dst[6] , src[6] );
00370         GF_ADDMULC( dst[7] , src[7] );
00371 #endif
00372 #if (UNROLL > 8)
00373         GF_ADDMULC( dst[8] , src[8] );
00374         GF_ADDMULC( dst[9] , src[9] );
00375         GF_ADDMULC( dst[10] , src[10] );
00376         GF_ADDMULC( dst[11] , src[11] );
00377         GF_ADDMULC( dst[12] , src[12] );
00378         GF_ADDMULC( dst[13] , src[13] );
00379         GF_ADDMULC( dst[14] , src[14] );
00380         GF_ADDMULC( dst[15] , src[15] );
00381 #endif
00382     }
00383 #endif
00384     lim += UNROLL - 1 ;
00385     for (; dst < lim; dst++, src++ )            /* final components */
00386         GF_ADDMULC( *dst , *src );
00387 }
00388 
00389 /*
00390  * computes C = AB where A is n*k, B is k*m, C is n*m
00391  */
00392 static void
00393 matmul(gf *a, gf *b, gf *c, int n, int k, int m)
00394 {
00395     int row, col, i ;
00396 
00397     for (row = 0; row < n ; row++) {
00398         for (col = 0; col < m ; col++) {
00399             gf *pa = &a[ row * k ];
00400             gf *pb = &b[ col ];
00401             gf acc = 0 ;
00402             for (i = 0; i < k ; i++, pa++, pb += m )
00403                 acc ^= gf_mul( *pa, *pb ) ;
00404             c[ row * m + col ] = acc ;
00405         }
00406     }
00407 }
00408 #ifdef DEBUG
00409 /*
00410  * returns 1 if the square matrix is identiy
00411  * (only for test)
00412  */
00413 /*
00414 static int
00415 is_identity(gf *m, int k)
00416 {
00417     int row, col ;
00418     for (row=0; row<k; row++)
00419         for (col=0; col<k; col++)
00420             if ( (row==col && *m != 1) ||
00421                  (row!=col && *m != 0) )
00422                  return 0 ;
00423             else
00424                 m++ ;
00425     return 1 ;
00426 }
00427 */
00428 #endif /* debug */
00429 
00430 /*
00431  * invert_mat() takes a matrix and produces its inverse
00432  * k is the size of the matrix.
00433  * (Gauss-Jordan, adapted from Numerical Recipes in C)
00434  * Return non-zero if singular.
00435  */
00436 DEB( int pivloops=0; int pivswaps=0 ; /* diagnostic */)
00437 static int
00438 invert_mat(gf *src, int k)
00439 {
00440     gf c, *p ;
00441     int irow, icol, row, col, i, ix ;
00442 
00443     int error = 1 ;
00444     int *indxc = (int*)my_malloc(k*sizeof(int), "indxc");
00445     int *indxr = (int*)my_malloc(k*sizeof(int), "indxr");
00446     int *ipiv = (int*)my_malloc(k*sizeof(int), "ipiv");
00447     gf *id_row = NEW_GF_MATRIX(1, k);
00448     gf *temp_row = NEW_GF_MATRIX(1, k);
00449 
00450     bzero(id_row, k*sizeof(gf));
00451     DEB( pivloops=0; pivswaps=0 ; /* diagnostic */ )
00452     /*
00453      * ipiv marks elements already used as pivots.
00454      */
00455     for (i = 0; i < k ; i++)
00456         ipiv[i] = 0 ;
00457 
00458     for (col = 0; col < k ; col++) {
00459         gf *pivot_row ;
00460         /*
00461          * Zeroing column 'col', look for a non-zero element.
00462          * First try on the diagonal, if it fails, look elsewhere.
00463          */
00464         irow = icol = -1 ;
00465         if (ipiv[col] != 1 && src[col*k + col] != 0) {
00466             irow = col ;
00467             icol = col ;
00468             goto found_piv ;
00469         }
00470         for (row = 0 ; row < k ; row++) {
00471             if (ipiv[row] != 1) {
00472                 for (ix = 0 ; ix < k ; ix++) {
00473                     DEB( pivloops++ ; )
00474                     if (ipiv[ix] == 0) {
00475                         if (src[row*k + ix] != 0) {
00476                             irow = row ;
00477                             icol = ix ;
00478                             goto found_piv ;
00479                         }
00480                     } else if (ipiv[ix] > 1) {
00481                         fprintf( stderr, "singular matrix\n" );
00482                         goto fail ;
00483                     }
00484                 }
00485             }
00486         }
00487         if (icol == -1) {
00488             fprintf( stderr, "XXX pivot not found!\n" );
00489             goto fail ;
00490         }
00491 found_piv:
00492         ++(ipiv[icol]) ;
00493         /*
00494          * swap rows irow and icol, so afterwards the diagonal
00495          * element will be correct. Rarely done, not worth
00496          * optimizing.
00497          */
00498         if (irow != icol) {
00499             for (ix = 0 ; ix < k ; ix++ ) {
00500                 SWAP( src[irow*k + ix], src[icol*k + ix], gf) ;
00501             }
00502         }
00503         indxr[col] = irow ;
00504         indxc[col] = icol ;
00505         pivot_row = &src[icol*k] ;
00506         c = pivot_row[icol] ;
00507         if (c == 0) {
00508             fprintf( stderr, "singular matrix 2\n" );
00509             goto fail ;
00510         }
00511         if (c != 1 ) { /* otherwhise this is a NOP */
00512             /*
00513              * this is done often , but optimizing is not so
00514              * fruitful, at least in the obvious ways (unrolling)
00515              */
00516             DEB( pivswaps++ ; )
00517             c = inverse[ c ] ;
00518             pivot_row[icol] = 1 ;
00519             for (ix = 0 ; ix < k ; ix++ )
00520                 pivot_row[ix] = gf_mul(c, pivot_row[ix] );
00521         }
00522         /*
00523          * from all rows, remove multiples of the selected row
00524          * to zero the relevant entry (in fact, the entry is not zero
00525          * because we know it must be zero).
00526          * (Here, if we know that the pivot_row is the identity,
00527          * we can optimize the addmul).
00528          */
00529         id_row[icol] = 1;
00530         if (bcmp(pivot_row, id_row, k*sizeof(gf)) != 0) {
00531             for (p = src, ix = 0 ; ix < k ; ix++, p += k ) {
00532                 if (ix != icol) {
00533                     c = p[icol] ;
00534                     p[icol] = 0 ;
00535                     addmul(p, pivot_row, c, k );
00536                 }
00537             }
00538         }
00539         id_row[icol] = 0;
00540     } /* done all columns */
00541     for (col = k-1 ; col >= 0 ; col-- ) {
00542         if (indxr[col] <0 || indxr[col] >= k)
00543            fprintf( stderr, "AARGH, indxr[col] %d\n", indxr[col] );
00544         else if (indxc[col] <0 || indxc[col] >= k)
00545             fprintf( stderr, "AARGH, indxc[col] %d\n", indxc[col] );
00546         else
00547         if (indxr[col] != indxc[col] ) {
00548             for (row = 0 ; row < k ; row++ ) {
00549                 SWAP( src[row*k + indxr[col]], src[row*k + indxc[col]], gf) ;
00550             }
00551         }
00552     }
00553     error = 0 ;
00554 fail:
00555     free(indxc);
00556     free(indxr);
00557     free(ipiv);
00558     free(id_row);
00559     free(temp_row);
00560     return error ;
00561 }
00562 
00563 /*
00564  * fast code for inverting a vandermonde matrix.
00565  * XXX NOTE: It assumes that the matrix
00566  * is not singular and _IS_ a vandermonde matrix. Only uses
00567  * the second column of the matrix, containing the p_i's.
00568  *
00569  * Algorithm borrowed from "Numerical recipes in C" -- sec.2.8, but
00570  * largely revised for my purposes.
00571  * p = coefficients of the matrix (p_i)
00572  * q = values of the polynomial (known)
00573  */
00574 
00575 int
00576 invert_vdm(gf *src, int k)
00577 {
00578     int i, j, row, col ;
00579     gf *b, *c, *p;
00580     gf t, xx ;
00581 
00582     if (k == 1)         /* degenerate case, matrix must be p^0 = 1 */
00583         return 0 ;
00584     /*
00585      * c holds the coefficient of P(x) = Prod (x - p_i), i=0..k-1
00586      * b holds the coefficient for the matrix inversion
00587      */
00588     c = NEW_GF_MATRIX(1, k);
00589     b = NEW_GF_MATRIX(1, k);
00590 
00591     p = NEW_GF_MATRIX(1, k);
00592 
00593     for ( j=1, i = 0 ; i < k ; i++, j+=k ) {
00594         c[i] = 0 ;
00595         p[i] = src[j] ;    /* p[i] */
00596     }
00597    /*
00598      * construct coeffs. recursively. We know c[k] = 1 (implicit)
00599      * and start P_0 = x - p_0, then at each stage multiply by
00600      * x - p_i generating P_i = x P_{i-1} - p_i P_{i-1}
00601      * After k steps we are done.
00602      */
00603     c[k-1] = p[0] ;     /* really -p(0), but x = -x in GF(2^m) */
00604     for (i = 1 ; i < k ; i++ ) {
00605         gf p_i = p[i] ; /* see above comment */
00606         for (j = k-1  - ( i - 1 ) ; j < k-1 ; j++ )
00607             c[j] ^= gf_mul( p_i, c[j+1] ) ;
00608         c[k-1] ^= p_i ;
00609     }
00610 
00611     for (row = 0 ; row < k ; row++ ) {
00612         /*
00613          * synthetic division etc.
00614          */
00615         xx = p[row] ;
00616         t = 1 ;
00617         b[k-1] = 1 ; /* this is in fact c[k] */
00618         for (i = k-2 ; i >= 0 ; i-- ) {
00619             b[i] = c[i+1] ^ gf_mul(xx, b[i+1]) ;
00620             t = gf_mul(xx, t) ^ b[i] ;
00621         }
00622         for (col = 0 ; col < k ; col++ )
00623             src[col*k + row] = gf_mul(inverse[t], b[col] );
00624     }
00625     free(c) ;
00626     free(b) ;
00627     free(p) ;
00628     return 0 ;
00629 }
00630 
00631 static int fec_initialized = 0 ;
00632 /* static */ void               /* VR: removed static */
00633 init_fec()
00634 {
00635     TICK(ticks[0]);
00636     generate_gf();
00637     TOCK(ticks[0]);
00638     DDB( fprintf( stderr, "generate_gf took %ldus\n", ticks[0] ); )
00639     TICK(ticks[0]);
00640     init_mul_table();
00641     TOCK(ticks[0]);
00642     DDB( fprintf( stderr, "init_mul_table took %ldus\n", ticks[0] ); )
00643     fec_initialized = 1 ;
00644 }
00645 
00646 /*
00647  * This section contains the proper FEC encoding/decoding routines.
00648  * The encoding matrix is computed starting with a Vandermonde matrix,
00649  * and then transforming it into a systematic matrix.
00650  */
00651 
00652 #define FEC_MAGIC       0xFECC0DEC
00653 
00654 struct fec_parms {
00655     u_long magic ;
00656     int k, n ;          /* parameters of the code */
00657     gf *enc_matrix ;
00658 } ;
00659 
00660 
00661 #define CPLUSPLUS_COMPATIBLE                            /* VR: added */
00662 #ifdef CPLUSPLUS_COMPATIBLE
00663 void fec_free(void *p_vp)
00664 #else
00665 void fec_free(struct fec_parms *p)
00666 #endif /* CPLUSPLUS_COMPATIBLE */
00667 {
00668 #ifdef CPLUSPLUS_COMPATIBLE
00669     struct fec_parms *p = (struct fec_parms *)p_vp;     /* VR */
00670 #endif /* CPLUSPLUS_COMPATIBLE */
00671     if (p==NULL ||
00672        p->magic != ( ( (FEC_MAGIC ^ p->k) ^ p->n) ^ (int)(p->enc_matrix)) ) {
00673         fprintf( stderr, "bad parameters to fec_free\n" );
00674         return ;
00675     }
00676     free(p->enc_matrix);
00677     free(p);
00678 }
00679 
00680 /*
00681  * create a new encoder, returning a descriptor. This contains k,n and
00682  * the encoding matrix.
00683  */
00684 #if 0           /* VR: changed as it creates problems with C++ compilers */
00685 struct fec_parms *
00686 #else
00687 void *
00688 #endif
00689 fec_new(int k, int n)
00690 {
00691     int row, col ;
00692     gf *p, *tmp_m ;
00693 
00694     struct fec_parms *retval ;
00695 
00696     if (fec_initialized == 0)
00697         init_fec();
00698 
00699     if (k > GF_SIZE + 1 || n > GF_SIZE + 1 || k > n ) {
00700         fprintf( stderr, "Invalid parameters k %d n %d GF_SIZE %d\n",
00701                 k, n, GF_SIZE );
00702         return NULL ;
00703     }
00704     retval = (struct fec_parms*)my_malloc(sizeof(struct fec_parms), "new_code");
00705     retval->k = k ;
00706     retval->n = n ;
00707     retval->enc_matrix = NEW_GF_MATRIX(n, k);
00708     retval->magic = ( ( FEC_MAGIC ^ k) ^ n) ^ (int)(retval->enc_matrix) ;
00709     tmp_m = NEW_GF_MATRIX(n, k);
00710     /*
00711      * fill the matrix with powers of field elements, starting from 0.
00712      * The first row is special, cannot be computed with exp. table.
00713      */
00714     tmp_m[0] = 1 ;
00715     for (col = 1; col < k ; col++)
00716         tmp_m[col] = 0 ;
00717     for (p = tmp_m + k, row = 0; row < n-1 ; row++, p += k) {
00718         for ( col = 0 ; col < k ; col ++ )
00719             p[col] = gf_exp[modnn(row*col)];
00720     }
00721 
00722     /*
00723      * quick code to build systematic matrix: invert the top
00724      * k*k vandermonde matrix, multiply right the bottom n-k rows
00725      * by the inverse, and construct the identity matrix at the top.
00726      */
00727     TICK(ticks[3]);
00728     invert_vdm(tmp_m, k); /* much faster than invert_mat */
00729     matmul(tmp_m + k*k, tmp_m, retval->enc_matrix + k*k, n - k, k, k);
00730     /*
00731      * the upper matrix is I so do not bother with a slow multiply
00732      */
00733     bzero(retval->enc_matrix, k*k*sizeof(gf) );
00734     for (p = retval->enc_matrix, col = 0 ; col < k ; col++, p += k+1 )
00735         *p = 1 ;
00736     free(tmp_m);
00737     TOCK(ticks[3]);
00738 
00739     DDB( fprintf( stderr, "--- %ld us to build encoding matrix\n", ticks[3] ); )
00740     DEB(pr_matrix(retval->enc_matrix, n, k, "encoding_matrix");)
00741 #if 0           /* VR: changed as it creates problems with C++ compilers */
00742     return retval ;
00743 #else
00744     return (void*)retval ;
00745 #endif
00746 }
00747 
00748 /*
00749  * fec_encode accepts as input pointers to n data packets of size sz,
00750  * and produces as output a packet pointed to by fec, computed
00751  * with index "index".
00752  */
00753 /*
00754  * VR: changed for C++ compilers who don't accept diff in parameters...
00755  * Use a definition that matches prototype in fec.h
00756  */
00757 #define CPLUSPLUS_COMPATIBLE                    /* VR: added */
00758 #ifdef CPLUSPLUS_COMPATIBLE
00759 void
00760 fec_encode(void *code_vp, void **src_vp, void *fec_vp, int index, int sz)
00761 #else
00762 void
00763 fec_encode(struct fec_parms *code, gf *src[], gf *fec, int index, int sz)
00764 #endif
00765 {
00766 #ifdef CPLUSPLUS_COMPATIBLE
00767     struct fec_parms *code = (struct fec_parms*)code_vp;/* VR */
00768     gf **src = (gf**)src_vp;                            /* VR */
00769     gf *fec = (gf*)fec_vp;                              /* VR */
00770 #endif /* CPLUSPLUS_COMPATIBLE */
00771     int i, k = code->k ;
00772     gf *p ;
00773 
00774     if (GF_BITS > 8)
00775         sz /= 2 ;
00776 
00777     if (index < k)
00778          bcopy(src[index], fec, sz*sizeof(gf) ) ;
00779     else if (index < code->n) {
00780         p = &(code->enc_matrix[index*k] );
00781         bzero(fec, sz*sizeof(gf));
00782         for (i = 0; i < k ; i++)
00783             addmul(fec, src[i], p[i], sz ) ;
00784     } else
00785         fprintf( stderr, "Invalid index %d (max %d)\n", index, code->n - 1 );
00786 }
00787 
00788 /*
00789  * shuffle move src packets in their position
00790  */
00791 static int
00792 shuffle(gf *pkt[], int index[], int k)
00793 {
00794     int i;
00795 
00796     for ( i = 0 ; i < k ; ) {
00797         if (index[i] >= k || index[i] == i)
00798             i++ ;
00799         else {
00800             /*
00801              * put pkt in the right position (first check for conflicts).
00802              */
00803             int c = index[i] ;
00804 
00805             if (index[c] == c) {
00806                 DEB( fprintf( stderr, "\nshuffle, error at %d\n", i ); )
00807                 return 1 ;
00808             }
00809             SWAP(index[i], index[c], int) ;
00810             SWAP(pkt[i], pkt[c], gf *) ;
00811         }
00812     }
00813     DEB( /* just test that it works... */
00814     for ( i = 0 ; i < k ; i++ ) {
00815         if (index[i] < k && index[i] != i) {
00816             fprintf( stderr, "shuffle: after\n" );
00817             for (i=0; i<k ; i++) fprintf( stderr, "%3d ", index[i] );
00818             fprintf( stderr, "\n" );
00819             return 1 ;
00820         }
00821     }
00822     )
00823     return 0 ;
00824 }
00825 
00826 /*
00827  * build_decode_matrix constructs the encoding matrix given the
00828  * indexes. The matrix must be already allocated as
00829  * a vector of k*k elements, in row-major order
00830  */
00831 static gf *
00832 build_decode_matrix(struct fec_parms *code, gf *pkt[], int index[])
00833 {
00834     int i , k = code->k ;
00835     gf *p, *matrix = NEW_GF_MATRIX(k, k);
00836 
00837     TICK(ticks[9]);
00838     for (i = 0, p = matrix ; i < k ; i++, p += k ) {
00839 #if 1 /* this is simply an optimization, not very useful indeed */
00840         if (index[i] < k) {
00841             bzero(p, k*sizeof(gf) );
00842             p[i] = 1 ;
00843         } else
00844 #endif
00845         if (index[i] < code->n )
00846             bcopy( &(code->enc_matrix[index[i]*k]), p, k*sizeof(gf) );
00847         else {
00848             fprintf( stderr, "decode: invalid index %d (max %d)\n",
00849                     index[i], code->n - 1 );
00850             free(matrix);
00851             return NULL ;
00852         }
00853     }
00854     TICK(ticks[9]);
00855     if (invert_mat(matrix, k)) {
00856         free(matrix);
00857         matrix = NULL ;
00858     }
00859     TOCK(ticks[9]);
00860     return matrix ;
00861 }
00862 
00863 /*
00864  * fec_decode receives as input a vector of packets, the indexes of
00865  * packets, and produces the correct vector as output.
00866  *
00867  *
00868  * Input:
00869  *      code: pointer to code descriptor
00870  *      pkt:  pointers to received packets. They are modified
00871  *            to store the output packets (in place)
00872  *      index: pointer to packet indexes (modified)
00873  *      sz:    size of each packet
00874  */
00875 /*
00876  * VR: changed for C++ compilers who don't accept diff in parameters...
00877  * Use a definition that matches prototype in fec.h
00878  */
00879 #define CPLUSPLUS_COMPATIBLE                    /* VR: added */
00880 #ifdef CPLUSPLUS_COMPATIBLE
00881 int
00882 fec_decode(void *code_vp, void **pkt_vp, int index[], int sz)
00883 #else
00884 int
00885 fec_decode(struct fec_parms *code, gf *pkt[], int index[], int sz)
00886 #endif
00887 {
00888 #ifdef CPLUSPLUS_COMPATIBLE
00889     struct fec_parms *code = (struct fec_parms*)code_vp;/* VR */
00890     gf **pkt = (gf**)pkt_vp;                            /* VR */
00891 #endif /* CPLUSPLUS_COMPATIBLE */
00892     gf *m_dec ;
00893     gf **new_pkt ;
00894     int row, col , k = code->k ;
00895 
00896     if (GF_BITS > 8)
00897         sz /= 2 ;
00898 
00899     if (shuffle(pkt, index, k)) /* error if TRUE */
00900         return 1 ;
00901     m_dec = build_decode_matrix(code, pkt, index);
00902 
00903     if (m_dec == NULL)
00904         return 1 ; /* error */
00905     /*
00906      * do the actual decoding
00907      */
00908     new_pkt = (gf **)my_malloc (k * sizeof (gf * ), "new pkt pointers" );
00909     for (row = 0 ; row < k ; row++ ) {
00910         if (index[row] >= k) {
00911             new_pkt[row] = (gf *)my_malloc (sz * sizeof (gf), "new pkt buffer" );
00912             bzero(new_pkt[row], sz * sizeof(gf) ) ;
00913             for (col = 0 ; col < k ; col++ )
00914                 addmul(new_pkt[row], pkt[col], m_dec[row*k + col], sz) ;
00915         }
00916     }
00917     /*
00918      * move pkts to their final destination
00919      */
00920     for (row = 0 ; row < k ; row++ ) {
00921         if (index[row] >= k) {
00922             bcopy(new_pkt[row], pkt[row], sz*sizeof(gf));
00923             free(new_pkt[row]);
00924         }
00925     }
00926     free(new_pkt);
00927     free(m_dec);
00928 
00929     return 0;
00930 }
00931 
00932 /*********** end of FEC code -- beginning of test code ************/
00933 
00934 #if (TEST || DEBUG)
00935 void
00936 test_gf()
00937 {
00938     int i ;
00939     /*
00940      * test gf tables. Sufficiently tested...
00941      */
00942     for (i=0; i<= GF_SIZE; i++) {
00943         if (gf_exp[gf_log[i]] != i)
00944             fprintf( stderr, "bad exp/log i %d log %d exp(log) %d\n",
00945                 i, gf_log[i], gf_exp[gf_log[i]]);
00946 
00947         if (i != 0 && gf_mul(i, inverse[i]) != 1)
00948             fprintf( stderr, "bad mul/inv i %d inv %d i*inv(i) %d\n",
00949                 i, inverse[i], gf_mul(i, inverse[i]) );
00950         if (gf_mul(0,i) != 0)
00951             fprintf( stderr, "bad mul table 0,%d\n",i);
00952         if (gf_mul(i,0) != 0)
00953             fprintf( stderr, "bad mul table %d,0\n",i);
00954     }
00955 }
00956 #endif /* TEST */
00957 
00958 #endif /* FEC */

Generated on Fri Mar 9 20:08:51 2007 for MAD-FCL by  doxygen 1.5.0