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