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