1 /** \file fec.c \brief Forward error correction based on Vandermonde matrices
2 * <br>980624<br> (C) 1997-98 Luigi Rizzo (luigi@iet.unipi.it)
3 *
4 * $Author: peltotal $ $Date: 2007/01/12 11:30:17 $ $Revision: 1.13 $
5 *
6 * Portions derived from code by Phil Karn (karn@ka9q.ampr.org),
7 * Robert Morelos-Zaragoza (robert@spectra.eng.hawaii.edu) and Hari
8 * Thirumoorthy (harit@spectra.eng.hawaii.edu), Aug 1995
9 *
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
12 * are met:
13 *
14 * 1. Redistributions of source code must retain the above copyright
15 * notice, this list of conditions and the following disclaimer.
16 * 2. Redistributions in binary form must reproduce the above
17 * copyright notice, this list of conditions and the following
18 * disclaimer in the documentation and/or other materials
19 * provided with the distribution.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND
22 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
23 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
24 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS
25 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
26 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
28 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
30 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
31 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
32 * OF SUCH DAMAGE.
33 */
34
35 #define FEC /* VR: added */
36
37 #ifdef FEC
38
39 #include "defines.h"
40 #include "fec.h"
41
42 #include <stdio.h>
43 #include <stdlib.h>
44 #include <string.h>
45
46 #ifndef _MSC_VER
47 #include <strings.h> /* VR: added */
48 #endif
49
50 #include <sys/types.h> /* VR: added */
51
52 /*
53 * compatibility stuff
54 */
55 #if (defined(MSDOS) || defined(_MSC_VER)) /* but also for others, e.g. sun... */
56 #define NEED_BCOPY
57 #define bcmp(a,b,n) memcmp(a,b,n)
58 #endif
59
60 #ifdef NEED_BCOPY
61 #define bcopy(s, d, siz) memcpy((d), (s), (siz))
62 #define bzero(d, siz) memset((d), '\0', (siz))
63 #endif
64
65 #ifndef u_long
66 #define u_long unsigned long
67 #endif
68
69 /*
70 * stuff used for testing purposes only
71 */
72
73 #ifdef TICK /* VR: avoid a warning under Solaris */
74 #undef TICK
75 #endif
76
77 #ifdef TEST
78 #define DEB(x)
79 #define DDB(x) x
80 #define DEBUG 4 /* minimal debugging */
81 #if (defined(MSDOS) || defined(_MSC_VER))
82 #include <time.h>
83 struct timeval {
84 unsigned long ticks;
85 };
86 #define gettimeofday(x, dummy) { (x)->ticks = clock() ; }
87 #define DIFF_T(a,b) (1+ 1000000*(a.ticks - b.ticks) / CLOCKS_PER_SEC )
88 typedef unsigned long u_long ;
89 typedef unsigned short u_short ;
90 #else /* typically, unix systems */
91 #include <sys/time.h>
92 #define DIFF_T(a,b) \
93 (1+ 1000000*(a.tv_sec - b.tv_sec) + (a.tv_usec - b.tv_usec) )
94 #endif
95
96 #define TICK(t) \
97 {struct timeval x ; \
98 gettimeofday(&x, NULL) ; \
99 t = x.tv_usec + 1000000* (x.tv_sec & 0xff ) ; \
100 }
101 #define TOCK(t) \
102 { u_long t1 ; TICK(t1) ; \
103 if (t1 < t) t = 256000000 + t1 - t ; \
104 else t = t1 - t ; \
105 if (t == 0) t = 1 ;}
106
107 u_long ticks[10]; /* vars for timekeeping */
108 #else
109 #define DEB(x)
110 #define DDB(x)
111 #define TICK(x)
112 #define TOCK(x)
113 #endif /* TEST */
114
115 /*
116 * You should not need to change anything beyond this point.
117 * The first part of the file implements linear algebra in GF.
118 *
119 * gf is the type used to store an element of the Galois Field.
120 * Must constain at least GF_BITS bits.
121 *
122 * Note: unsigned char will work up to GF(256) but int seems to run
123 * faster on the Pentium. We use int whenever have to deal with an
124 * index, since they are generally faster.
125 */
126 #if (GF_BITS < 2 || GF_BITS >16)
127 #error "GF_BITS must be 2 .. 16"
128 #endif
129
130 #if (GF_BITS <= 8)
131 typedef unsigned char gf;
132 #else
133 typedef unsigned short gf;
134 #endif
135
136 // defined in fec.h
137 //#define GF_SIZE ((1 << GF_BITS) - 1) /* powers of \alpha */
138
139 /*
140 * Primitive polynomials - see Lin & Costello, Appendix A,
141 * and Lee & Messerschmitt, p. 453.
142 */
143 static char *allPp[] = { /* GF_BITS polynomial */
144 NULL, /* 0 no code */
145 NULL, /* 1 no code */
146 "111", /* 2 1+x+x^2 */
147 "1101", /* 3 1+x+x^3 */
148 "11001", /* 4 1+x+x^4 */
149 "101001", /* 5 1+x^2+x^5 */
150 "1100001", /* 6 1+x+x^6 */
151 "10010001", /* 7 1 + x^3 + x^7 */
152 "101110001", /* 8 1+x^2+x^3+x^4+x^8 */
153 "1000100001", /* 9 1+x^4+x^9 */
154 "10010000001", /* 10 1+x^3+x^10 */
155 "101000000001", /* 11 1+x^2+x^11 */
156 "1100101000001", /* 12 1+x+x^4+x^6+x^12 */
157 "11011000000001", /* 13 1+x+x^3+x^4+x^13 */
158 "110000100010001", /* 14 1+x+x^6+x^10+x^14 */
159 "1100000000000001", /* 15 1+x+x^15 */
160 "11010000000010001" /* 16 1+x+x^3+x^12+x^16 */
161 };
162
163
164 /*
165 * To speed up computations, we have tables for logarithm, exponent
166 * and inverse of a number. If GF_BITS <= 8, we use a table for
167 * multiplication as well (it takes 64K, no big deal even on a PDA,
168 * especially because it can be pre-initialized an put into a ROM!),
169 * otherwhise we use a table of logarithms.
170 * In any case the macro gf_mul(x,y) takes care of multiplications.
171 */
172
173 static gf gf_exp[2*GF_SIZE]; /* index->poly form conversion table */
174 static int gf_log[GF_SIZE + 1]; /* Poly->index form conversion table */
175 static gf inverse[GF_SIZE+1]; /* inverse of field elem. */
176 /* inv[\alpha**i]=\alpha**(GF_SIZE-i-1) */
177
178 /*
179 * modnn(x) computes x % GF_SIZE, where GF_SIZE is 2**GF_BITS - 1,
180 * without a slow divide.
181 */
182 static gf
modnn(int x)183 modnn(int x)
184 {
185 while (x >= GF_SIZE) {
186 x -= GF_SIZE;
187 x = (x >> GF_BITS) + (x & GF_SIZE);
188 }
189 return x;
190 }
191
192 #define SWAP(a,b,t) {t tmp; tmp=a; a=b; b=tmp;}
193
194 /*
195 * gf_mul(x,y) multiplies two numbers. If GF_BITS<=8, it is much
196 * faster to use a multiplication table.
197 *
198 * USE_GF_MULC, GF_MULC0(c) and GF_ADDMULC(x) can be used when multiplying
199 * many numbers by the same constant. In this case the first
200 * call sets the constant, and others perform the multiplications.
201 * A value related to the multiplication is held in a local variable
202 * declared with USE_GF_MULC . See usage in addmul1().
203 */
204 #if (GF_BITS <= 8)
205 static gf gf_mul_table[GF_SIZE + 1][GF_SIZE + 1];
206
207 #define gf_mul(x,y) gf_mul_table[x][y]
208
209 #define USE_GF_MULC register gf * __gf_mulc_
210 #define GF_MULC0(c) __gf_mulc_ = gf_mul_table[c]
211 #define GF_ADDMULC(dst, x) dst ^= __gf_mulc_[x]
212
213 static void
init_mul_table()214 init_mul_table()
215 {
216 int i, j;
217 for (i=0; i< GF_SIZE+1; i++)
218 for (j=0; j< GF_SIZE+1; j++)
219 gf_mul_table[i][j] = gf_exp[modnn(gf_log[i] + gf_log[j]) ] ;
220
221 for (j=0; j< GF_SIZE+1; j++)
222 gf_mul_table[0][j] = gf_mul_table[j][0] = 0;
223 }
224 #else /* GF_BITS > 8 */
225 static inline gf
gf_mul(x,y)226 gf_mul(x,y)
227 {
228 if ( (x) == 0 || (y)==0 ) return 0;
229
230 return gf_exp[gf_log[x] + gf_log[y] ] ;
231 }
232 #define init_mul_table()
233
234 #define USE_GF_MULC register gf * __gf_mulc_
235 #define GF_MULC0(c) __gf_mulc_ = &gf_exp[ gf_log[c] ]
236 #define GF_ADDMULC(dst, x) { if (x) dst ^= __gf_mulc_[ gf_log[x] ] ; }
237 #endif
238
239 /*
240 * Generate GF(2**m) from the irreducible polynomial p(X) in p[0]..p[m]
241 * Lookup tables:
242 * index->polynomial form gf_exp[] contains j= \alpha^i;
243 * polynomial form -> index form gf_log[ j = \alpha^i ] = i
244 * \alpha=x is the primitive element of GF(2^m)
245 *
246 * For efficiency, gf_exp[] has size 2*GF_SIZE, so that a simple
247 * multiplication of two numbers can be resolved without calling modnn
248 */
249
250 /*
251 * i use malloc so many times, it is easier to put checks all in
252 * one place.
253 */
254 static void *
my_malloc(int sz,char * err_string)255 my_malloc(int sz, char *err_string)
256 {
257 void *p = malloc( sz );
258 if (p == NULL) {
259 fprintf( stderr, "-- malloc failure allocation %s\n", err_string );
260 exit(1) ;
261 }
262 return p ;
263 }
264
265 #define NEW_GF_MATRIX(rows, cols) \
266 (gf *)my_malloc(rows * cols * sizeof(gf), " ## __LINE__ ## " )
267
268 /*
269 * initialize the data structures used for computations in GF.
270 */
271 static void
generate_gf(void)272 generate_gf(void)
273 {
274 int i;
275 gf mask;
276 char *Pp = allPp[GF_BITS] ;
277
278 mask = 1; /* x ** 0 = 1 */
279 gf_exp[GF_BITS] = 0; /* will be updated at the end of the 1st loop */
280 /*
281 * first, generate the (polynomial representation of) powers of \alpha,
282 * which are stored in gf_exp[i] = \alpha ** i .
283 * At the same time build gf_log[gf_exp[i]] = i .
284 * The first GF_BITS powers are simply bits shifted to the left.
285 */
286 for (i = 0; i < GF_BITS; i++, mask <<= 1 ) {
287 gf_exp[i] = mask;
288 gf_log[gf_exp[i]] = i;
289 /*
290 * If Pp[i] == 1 then \alpha ** i occurs in poly-repr
291 * gf_exp[GF_BITS] = \alpha ** GF_BITS
292 */
293 if ( Pp[i] == '1' )
294 gf_exp[GF_BITS] ^= mask;
295 }
296 /*
297 * now gf_exp[GF_BITS] = \alpha ** GF_BITS is complete, so can als
298 * compute its inverse.
299 */
300 gf_log[gf_exp[GF_BITS]] = GF_BITS;
301 /*
302 * Poly-repr of \alpha ** (i+1) is given by poly-repr of
303 * \alpha ** i shifted left one-bit and accounting for any
304 * \alpha ** GF_BITS term that may occur when poly-repr of
305 * \alpha ** i is shifted.
306 */
307 mask = 1 << (GF_BITS - 1 ) ;
308 for (i = GF_BITS + 1; i < GF_SIZE; i++) {
309 if (gf_exp[i - 1] >= mask)
310 gf_exp[i] = gf_exp[GF_BITS] ^ ((gf_exp[i - 1] ^ mask) << 1);
311 else
312 gf_exp[i] = gf_exp[i - 1] << 1;
313 gf_log[gf_exp[i]] = i;
314 }
315 /*
316 * log(0) is not defined, so use a special value
317 */
318 gf_log[0] = GF_SIZE ;
319 /* set the extended gf_exp values for fast multiply */
320 for (i = 0 ; i < GF_SIZE ; i++)
321 gf_exp[i + GF_SIZE] = gf_exp[i] ;
322
323 /*
324 * again special cases. 0 has no inverse. This used to
325 * be initialized to GF_SIZE, but it should make no difference
326 * since noone is supposed to read from here.
327 */
328 inverse[0] = 0 ;
329 inverse[1] = 1;
330 for (i=2; i<=GF_SIZE; i++)
331 inverse[i] = gf_exp[GF_SIZE-gf_log[i]];
332 }
333
334 /*
335 * Various linear algebra operations that i use often.
336 */
337
338 /*
339 * addmul() computes dst[] = dst[] + c * src[]
340 * This is used often, so better optimize it! Currently the loop is
341 * unrolled 16 times, a good value for 486 and pentium-class machines.
342 * The case c=0 is also optimized, whereas c=1 is not. These
343 * calls are unfrequent in my typical apps so I did not bother.
344 *
345 * Note that gcc on
346 */
347 #define addmul(dst, src, c, sz) \
348 if (c != 0) addmul1(dst, src, c, sz)
349
350 #define UNROLL 16 /* 1, 4, 8, 16 */
351 static void
addmul1(gf * dst1,gf * src1,gf c,int sz)352 addmul1(gf *dst1, gf *src1, gf c, int sz)
353 {
354 USE_GF_MULC ;
355 register gf *dst = dst1, *src = src1 ;
356 gf *lim = &dst[sz - UNROLL + 1] ;
357
358 GF_MULC0(c) ;
359
360 #if (UNROLL > 1) /* unrolling by 8/16 is quite effective on the pentium */
361 for (; dst < lim ; dst += UNROLL, src += UNROLL ) {
362 GF_ADDMULC( dst[0] , src[0] );
363 GF_ADDMULC( dst[1] , src[1] );
364 GF_ADDMULC( dst[2] , src[2] );
365 GF_ADDMULC( dst[3] , src[3] );
366 #if (UNROLL > 4)
367 GF_ADDMULC( dst[4] , src[4] );
368 GF_ADDMULC( dst[5] , src[5] );
369 GF_ADDMULC( dst[6] , src[6] );
370 GF_ADDMULC( dst[7] , src[7] );
371 #endif
372 #if (UNROLL > 8)
373 GF_ADDMULC( dst[8] , src[8] );
374 GF_ADDMULC( dst[9] , src[9] );
375 GF_ADDMULC( dst[10] , src[10] );
376 GF_ADDMULC( dst[11] , src[11] );
377 GF_ADDMULC( dst[12] , src[12] );
378 GF_ADDMULC( dst[13] , src[13] );
379 GF_ADDMULC( dst[14] , src[14] );
380 GF_ADDMULC( dst[15] , src[15] );
381 #endif
382 }
383 #endif
384 lim += UNROLL - 1 ;
385 for (; dst < lim; dst++, src++ ) /* final components */
386 GF_ADDMULC( *dst , *src );
387 }
388
389 /*
390 * computes C = AB where A is n*k, B is k*m, C is n*m
391 */
392 static void
matmul(gf * a,gf * b,gf * c,int n,int k,int m)393 matmul(gf *a, gf *b, gf *c, int n, int k, int m)
394 {
395 int row, col, i ;
396
397 for (row = 0; row < n ; row++) {
398 for (col = 0; col < m ; col++) {
399 gf *pa = &a[ row * k ];
400 gf *pb = &b[ col ];
401 gf acc = 0 ;
402 for (i = 0; i < k ; i++, pa++, pb += m )
403 acc ^= gf_mul( *pa, *pb ) ;
404 c[ row * m + col ] = acc ;
405 }
406 }
407 }
408 #ifdef DEBUG
409 /*
410 * returns 1 if the square matrix is identiy
411 * (only for test)
412 */
413 /*
414 static int
415 is_identity(gf *m, int k)
416 {
417 int row, col ;
418 for (row=0; row<k; row++)
419 for (col=0; col<k; col++)
420 if ( (row==col && *m != 1) ||
421 (row!=col && *m != 0) )
422 return 0 ;
423 else
424 m++ ;
425 return 1 ;
426 }
427 */
428 #endif /* debug */
429
430 /*
431 * invert_mat() takes a matrix and produces its inverse
432 * k is the size of the matrix.
433 * (Gauss-Jordan, adapted from Numerical Recipes in C)
434 * Return non-zero if singular.
435 */
436 DEB( int pivloops=0; int pivswaps=0 ; /* diagnostic */)
437 static int
invert_mat(gf * src,int k)438 invert_mat(gf *src, int k)
439 {
440 gf c, *p ;
441 int irow, icol, row, col, i, ix ;
442
443 int error = 1 ;
444 int *indxc = (int*)my_malloc(k*sizeof(int), "indxc");
445 int *indxr = (int*)my_malloc(k*sizeof(int), "indxr");
446 int *ipiv = (int*)my_malloc(k*sizeof(int), "ipiv");
447 gf *id_row = NEW_GF_MATRIX(1, k);
448 gf *temp_row = NEW_GF_MATRIX(1, k);
449
450 bzero(id_row, k*sizeof(gf));
451 DEB( pivloops=0; pivswaps=0 ; /* diagnostic */ )
452 /*
453 * ipiv marks elements already used as pivots.
454 */
455 for (i = 0; i < k ; i++)
456 ipiv[i] = 0 ;
457
458 for (col = 0; col < k ; col++) {
459 gf *pivot_row ;
460 /*
461 * Zeroing column 'col', look for a non-zero element.
462 * First try on the diagonal, if it fails, look elsewhere.
463 */
464 irow = icol = -1 ;
465 if (ipiv[col] != 1 && src[col*k + col] != 0) {
466 irow = col ;
467 icol = col ;
468 goto found_piv ;
469 }
470 for (row = 0 ; row < k ; row++) {
471 if (ipiv[row] != 1) {
472 for (ix = 0 ; ix < k ; ix++) {
473 DEB( pivloops++ ; )
474 if (ipiv[ix] == 0) {
475 if (src[row*k + ix] != 0) {
476 irow = row ;
477 icol = ix ;
478 goto found_piv ;
479 }
480 } else if (ipiv[ix] > 1) {
481 fprintf( stderr, "singular matrix\n" );
482 goto fail ;
483 }
484 }
485 }
486 }
487 if (icol == -1) {
488 fprintf( stderr, "XXX pivot not found!\n" );
489 goto fail ;
490 }
491 found_piv:
492 ++(ipiv[icol]) ;
493 /*
494 * swap rows irow and icol, so afterwards the diagonal
495 * element will be correct. Rarely done, not worth
496 * optimizing.
497 */
498 if (irow != icol) {
499 for (ix = 0 ; ix < k ; ix++ ) {
500 SWAP( src[irow*k + ix], src[icol*k + ix], gf) ;
501 }
502 }
503 indxr[col] = irow ;
504 indxc[col] = icol ;
505 pivot_row = &src[icol*k] ;
506 c = pivot_row[icol] ;
507 if (c == 0) {
508 fprintf( stderr, "singular matrix 2\n" );
509 goto fail ;
510 }
511 if (c != 1 ) { /* otherwhise this is a NOP */
512 /*
513 * this is done often , but optimizing is not so
514 * fruitful, at least in the obvious ways (unrolling)
515 */
516 DEB( pivswaps++ ; )
517 c = inverse[ c ] ;
518 pivot_row[icol] = 1 ;
519 for (ix = 0 ; ix < k ; ix++ )
520 pivot_row[ix] = gf_mul(c, pivot_row[ix] );
521 }
522 /*
523 * from all rows, remove multiples of the selected row
524 * to zero the relevant entry (in fact, the entry is not zero
525 * because we know it must be zero).
526 * (Here, if we know that the pivot_row is the identity,
527 * we can optimize the addmul).
528 */
529 id_row[icol] = 1;
530 if (bcmp(pivot_row, id_row, k*sizeof(gf)) != 0) {
531 for (p = src, ix = 0 ; ix < k ; ix++, p += k ) {
532 if (ix != icol) {
533 c = p[icol] ;
534 p[icol] = 0 ;
535 addmul(p, pivot_row, c, k );
536 }
537 }
538 }
539 id_row[icol] = 0;
540 } /* done all columns */
541 for (col = k-1 ; col >= 0 ; col-- ) {
542 if (indxr[col] <0 || indxr[col] >= k)
543 fprintf( stderr, "AARGH, indxr[col] %d\n", indxr[col] );
544 else if (indxc[col] <0 || indxc[col] >= k)
545 fprintf( stderr, "AARGH, indxc[col] %d\n", indxc[col] );
546 else
547 if (indxr[col] != indxc[col] ) {
548 for (row = 0 ; row < k ; row++ ) {
549 SWAP( src[row*k + indxr[col]], src[row*k + indxc[col]], gf) ;
550 }
551 }
552 }
553 error = 0 ;
554 fail:
555 free(indxc);
556 free(indxr);
557 free(ipiv);
558 free(id_row);
559 free(temp_row);
560 return error ;
561 }
562
563 /*
564 * fast code for inverting a vandermonde matrix.
565 * XXX NOTE: It assumes that the matrix
566 * is not singular and _IS_ a vandermonde matrix. Only uses
567 * the second column of the matrix, containing the p_i's.
568 *
569 * Algorithm borrowed from "Numerical recipes in C" -- sec.2.8, but
570 * largely revised for my purposes.
571 * p = coefficients of the matrix (p_i)
572 * q = values of the polynomial (known)
573 */
574
575 int
invert_vdm(gf * src,int k)576 invert_vdm(gf *src, int k)
577 {
578 int i, j, row, col ;
579 gf *b, *c, *p;
580 gf t, xx ;
581
582 if (k == 1) /* degenerate case, matrix must be p^0 = 1 */
583 return 0 ;
584 /*
585 * c holds the coefficient of P(x) = Prod (x - p_i), i=0..k-1
586 * b holds the coefficient for the matrix inversion
587 */
588 c = NEW_GF_MATRIX(1, k);
589 b = NEW_GF_MATRIX(1, k);
590
591 p = NEW_GF_MATRIX(1, k);
592
593 for ( j=1, i = 0 ; i < k ; i++, j+=k ) {
594 c[i] = 0 ;
595 p[i] = src[j] ; /* p[i] */
596 }
597 /*
598 * construct coeffs. recursively. We know c[k] = 1 (implicit)
599 * and start P_0 = x - p_0, then at each stage multiply by
600 * x - p_i generating P_i = x P_{i-1} - p_i P_{i-1}
601 * After k steps we are done.
602 */
603 c[k-1] = p[0] ; /* really -p(0), but x = -x in GF(2^m) */
604 for (i = 1 ; i < k ; i++ ) {
605 gf p_i = p[i] ; /* see above comment */
606 for (j = k-1 - ( i - 1 ) ; j < k-1 ; j++ )
607 c[j] ^= gf_mul( p_i, c[j+1] ) ;
608 c[k-1] ^= p_i ;
609 }
610
611 for (row = 0 ; row < k ; row++ ) {
612 /*
613 * synthetic division etc.
614 */
615 xx = p[row] ;
616 t = 1 ;
617 b[k-1] = 1 ; /* this is in fact c[k] */
618 for (i = k-2 ; i >= 0 ; i-- ) {
619 b[i] = c[i+1] ^ gf_mul(xx, b[i+1]) ;
620 t = gf_mul(xx, t) ^ b[i] ;
621 }
622 for (col = 0 ; col < k ; col++ )
623 src[col*k + row] = gf_mul(inverse[t], b[col] );
624 }
625 free(c) ;
626 free(b) ;
627 free(p) ;
628 return 0 ;
629 }
630
631 static int fec_initialized = 0 ;
632 /* static */ void /* VR: removed static */
init_fec()633 init_fec()
634 {
635 TICK(ticks[0]);
636 generate_gf();
637 TOCK(ticks[0]);
638 DDB( fprintf( stderr, "generate_gf took %ldus\n", ticks[0] ); )
639 TICK(ticks[0]);
640 init_mul_table();
641 TOCK(ticks[0]);
642 DDB( fprintf( stderr, "init_mul_table took %ldus\n", ticks[0] ); )
643 fec_initialized = 1 ;
644 }
645
646 /*
647 * This section contains the proper FEC encoding/decoding routines.
648 * The encoding matrix is computed starting with a Vandermonde matrix,
649 * and then transforming it into a systematic matrix.
650 */
651
652 #define FEC_MAGIC 0xFECC0DEC
653
654 struct fec_parms {
655 u_long magic ;
656 int k, n ; /* parameters of the code */
657 gf *enc_matrix ;
658 } ;
659
660
661 #define CPLUSPLUS_COMPATIBLE /* VR: added */
662 #ifdef CPLUSPLUS_COMPATIBLE
fec_free(void * p_vp)663 void fec_free(void *p_vp)
664 #else
665 void fec_free(struct fec_parms *p)
666 #endif /* CPLUSPLUS_COMPATIBLE */
667 {
668 #ifdef CPLUSPLUS_COMPATIBLE
669 struct fec_parms *p = (struct fec_parms *)p_vp; /* VR */
670 #endif /* CPLUSPLUS_COMPATIBLE */
671 if (p==NULL ||
672 p->magic != ( ( (FEC_MAGIC ^ p->k) ^ p->n) ^ (int)(p->enc_matrix)) ) {
673 fprintf( stderr, "bad parameters to fec_free\n" );
674 return ;
675 }
676 free(p->enc_matrix);
677 free(p);
678 }
679
680 /*
681 * create a new encoder, returning a descriptor. This contains k,n and
682 * the encoding matrix.
683 */
684 #if 0 /* VR: changed as it creates problems with C++ compilers */
685 struct fec_parms *
686 #else
687 void *
688 #endif
fec_new(int k,int n)689 fec_new(int k, int n)
690 {
691 int row, col ;
692 gf *p, *tmp_m ;
693
694 struct fec_parms *retval ;
695
696 if (fec_initialized == 0)
697 init_fec();
698
699 if (k > GF_SIZE + 1 || n > GF_SIZE + 1 || k > n ) {
700 fprintf( stderr, "Invalid parameters k %d n %d GF_SIZE %d\n",
701 k, n, GF_SIZE );
702 return NULL ;
703 }
704 retval = (struct fec_parms*)my_malloc(sizeof(struct fec_parms), "new_code");
705 retval->k = k ;
706 retval->n = n ;
707 retval->enc_matrix = NEW_GF_MATRIX(n, k);
708 retval->magic = ( ( FEC_MAGIC ^ k) ^ n) ^ (int)(retval->enc_matrix) ;
709 tmp_m = NEW_GF_MATRIX(n, k);
710 /*
711 * fill the matrix with powers of field elements, starting from 0.
712 * The first row is special, cannot be computed with exp. table.
713 */
714 tmp_m[0] = 1 ;
715 for (col = 1; col < k ; col++)
716 tmp_m[col] = 0 ;
717 for (p = tmp_m + k, row = 0; row < n-1 ; row++, p += k) {
718 for ( col = 0 ; col < k ; col ++ )
719 p[col] = gf_exp[modnn(row*col)];
720 }
721
722 /*
723 * quick code to build systematic matrix: invert the top
724 * k*k vandermonde matrix, multiply right the bottom n-k rows
725 * by the inverse, and construct the identity matrix at the top.
726 */
727 TICK(ticks[3]);
728 invert_vdm(tmp_m, k); /* much faster than invert_mat */
729 matmul(tmp_m + k*k, tmp_m, retval->enc_matrix + k*k, n - k, k, k);
730 /*
731 * the upper matrix is I so do not bother with a slow multiply
732 */
733 bzero(retval->enc_matrix, k*k*sizeof(gf) );
734 for (p = retval->enc_matrix, col = 0 ; col < k ; col++, p += k+1 )
735 *p = 1 ;
736 free(tmp_m);
737 TOCK(ticks[3]);
738
739 DDB( fprintf( stderr, "--- %ld us to build encoding matrix\n", ticks[3] ); )
740 DEB(pr_matrix(retval->enc_matrix, n, k, "encoding_matrix");)
741 #if 0 /* VR: changed as it creates problems with C++ compilers */
742 return retval ;
743 #else
744 return (void*)retval ;
745 #endif
746 }
747
748 /*
749 * fec_encode accepts as input pointers to n data packets of size sz,
750 * and produces as output a packet pointed to by fec, computed
751 * with index "index".
752 */
753 /*
754 * VR: changed for C++ compilers who don't accept diff in parameters...
755 * Use a definition that matches prototype in fec.h
756 */
757 #define CPLUSPLUS_COMPATIBLE /* VR: added */
758 #ifdef CPLUSPLUS_COMPATIBLE
759 void
fec_encode(void * code_vp,void ** src_vp,void * fec_vp,int index,int sz)760 fec_encode(void *code_vp, void **src_vp, void *fec_vp, int index, int sz)
761 #else
762 void
763 fec_encode(struct fec_parms *code, gf *src[], gf *fec, int index, int sz)
764 #endif
765 {
766 #ifdef CPLUSPLUS_COMPATIBLE
767 struct fec_parms *code = (struct fec_parms*)code_vp;/* VR */
768 gf **src = (gf**)src_vp; /* VR */
769 gf *fec = (gf*)fec_vp; /* VR */
770 #endif /* CPLUSPLUS_COMPATIBLE */
771 int i, k = code->k ;
772 gf *p ;
773
774 if (GF_BITS > 8)
775 sz /= 2 ;
776
777 if (index < k)
778 bcopy(src[index], fec, sz*sizeof(gf) ) ;
779 else if (index < code->n) {
780 p = &(code->enc_matrix[index*k] );
781 bzero(fec, sz*sizeof(gf));
782 for (i = 0; i < k ; i++)
783 addmul(fec, src[i], p[i], sz ) ;
784 } else
785 fprintf( stderr, "Invalid index %d (max %d)\n", index, code->n - 1 );
786 }
787
788 /*
789 * shuffle move src packets in their position
790 */
791 static int
shuffle(gf * pkt[],int index[],int k)792 shuffle(gf *pkt[], int index[], int k)
793 {
794 int i;
795
796 for ( i = 0 ; i < k ; ) {
797 if (index[i] >= k || index[i] == i)
798 i++ ;
799 else {
800 /*
801 * put pkt in the right position (first check for conflicts).
802 */
803 int c = index[i] ;
804
805 if (index[c] == c) {
806 DEB( fprintf( stderr, "\nshuffle, error at %d\n", i ); )
807 return 1 ;
808 }
809 SWAP(index[i], index[c], int) ;
810 SWAP(pkt[i], pkt[c], gf *) ;
811 }
812 }
813 DEB( /* just test that it works... */
814 for ( i = 0 ; i < k ; i++ ) {
815 if (index[i] < k && index[i] != i) {
816 fprintf( stderr, "shuffle: after\n" );
817 for (i=0; i<k ; i++) fprintf( stderr, "%3d ", index[i] );
818 fprintf( stderr, "\n" );
819 return 1 ;
820 }
821 }
822 )
823 return 0 ;
824 }
825
826 /*
827 * build_decode_matrix constructs the encoding matrix given the
828 * indexes. The matrix must be already allocated as
829 * a vector of k*k elements, in row-major order
830 */
831 static gf *
build_decode_matrix(struct fec_parms * code,gf * pkt[],int index[])832 build_decode_matrix(struct fec_parms *code, gf *pkt[], int index[])
833 {
834 int i , k = code->k ;
835 gf *p, *matrix = NEW_GF_MATRIX(k, k);
836
837 TICK(ticks[9]);
838 for (i = 0, p = matrix ; i < k ; i++, p += k ) {
839 #if 1 /* this is simply an optimization, not very useful indeed */
840 if (index[i] < k) {
841 bzero(p, k*sizeof(gf) );
842 p[i] = 1 ;
843 } else
844 #endif
845 if (index[i] < code->n )
846 bcopy( &(code->enc_matrix[index[i]*k]), p, k*sizeof(gf) );
847 else {
848 fprintf( stderr, "decode: invalid index %d (max %d)\n",
849 index[i], code->n - 1 );
850 free(matrix);
851 return NULL ;
852 }
853 }
854 TICK(ticks[9]);
855 if (invert_mat(matrix, k)) {
856 free(matrix);
857 matrix = NULL ;
858 }
859 TOCK(ticks[9]);
860 return matrix ;
861 }
862
863 /*
864 * fec_decode receives as input a vector of packets, the indexes of
865 * packets, and produces the correct vector as output.
866 *
867 *
868 * Input:
869 * code: pointer to code descriptor
870 * pkt: pointers to received packets. They are modified
871 * to store the output packets (in place)
872 * index: pointer to packet indexes (modified)
873 * sz: size of each packet
874 */
875 /*
876 * VR: changed for C++ compilers who don't accept diff in parameters...
877 * Use a definition that matches prototype in fec.h
878 */
879 #define CPLUSPLUS_COMPATIBLE /* VR: added */
880 #ifdef CPLUSPLUS_COMPATIBLE
881 int
fec_decode(void * code_vp,void ** pkt_vp,int index[],int sz)882 fec_decode(void *code_vp, void **pkt_vp, int index[], int sz)
883 #else
884 int
885 fec_decode(struct fec_parms *code, gf *pkt[], int index[], int sz)
886 #endif
887 {
888 #ifdef CPLUSPLUS_COMPATIBLE
889 struct fec_parms *code = (struct fec_parms*)code_vp;/* VR */
890 gf **pkt = (gf**)pkt_vp; /* VR */
891 #endif /* CPLUSPLUS_COMPATIBLE */
892 gf *m_dec ;
893 gf **new_pkt ;
894 int row, col , k = code->k ;
895
896 if (GF_BITS > 8)
897 sz /= 2 ;
898
899 if (shuffle(pkt, index, k)) /* error if TRUE */
900 return 1 ;
901 m_dec = build_decode_matrix(code, pkt, index);
902
903 if (m_dec == NULL)
904 return 1 ; /* error */
905 /*
906 * do the actual decoding
907 */
908 new_pkt = (gf **)my_malloc (k * sizeof (gf * ), "new pkt pointers" );
909 for (row = 0 ; row < k ; row++ ) {
910 if (index[row] >= k) {
911 new_pkt[row] = (gf *)my_malloc (sz * sizeof (gf), "new pkt buffer" );
912 bzero(new_pkt[row], sz * sizeof(gf) ) ;
913 for (col = 0 ; col < k ; col++ )
914 addmul(new_pkt[row], pkt[col], m_dec[row*k + col], sz) ;
915 }
916 }
917 /*
918 * move pkts to their final destination
919 */
920 for (row = 0 ; row < k ; row++ ) {
921 if (index[row] >= k) {
922 bcopy(new_pkt[row], pkt[row], sz*sizeof(gf));
923 free(new_pkt[row]);
924 }
925 }
926 free(new_pkt);
927 free(m_dec);
928
929 return 0;
930 }
931
932 /*********** end of FEC code -- beginning of test code ************/
933
934 #if (TEST || DEBUG)
935 void
test_gf()936 test_gf()
937 {
938 int i ;
939 /*
940 * test gf tables. Sufficiently tested...
941 */
942 for (i=0; i<= GF_SIZE; i++) {
943 if (gf_exp[gf_log[i]] != i)
944 fprintf( stderr, "bad exp/log i %d log %d exp(log) %d\n",
945 i, gf_log[i], gf_exp[gf_log[i]]);
946
947 if (i != 0 && gf_mul(i, inverse[i]) != 1)
948 fprintf( stderr, "bad mul/inv i %d inv %d i*inv(i) %d\n",
949 i, inverse[i], gf_mul(i, inverse[i]) );
950 if (gf_mul(0,i) != 0)
951 fprintf( stderr, "bad mul table 0,%d\n",i);
952 if (gf_mul(i,0) != 0)
953 fprintf( stderr, "bad mul table %d,0\n",i);
954 }
955 }
956 #endif /* TEST */
957
958 #endif /* FEC */
959