|
@@ -0,0 +1,998 @@
|
|
|
|
|
+/*#define PROFILE*/
|
|
|
|
|
+/*
|
|
|
|
|
+ * fec.c -- forward error correction based on Vandermonde matrices
|
|
|
|
|
+ * 980624
|
|
|
|
|
+ * (C) 1997-98 Luigi Rizzo (luigi@iet.unipi.it)
|
|
|
|
|
+ * (C) 2001 Alain Knaff (alain@knaff.lu)
|
|
|
|
|
+ *
|
|
|
|
|
+ * Portions derived from code by Phil Karn (karn@ka9q.ampr.org),
|
|
|
|
|
+ * Robert Morelos-Zaragoza (robert@spectra.eng.hawaii.edu) and Hari
|
|
|
|
|
+ * Thirumoorthy (harit@spectra.eng.hawaii.edu), Aug 1995
|
|
|
|
|
+ *
|
|
|
|
|
+ * Redistribution and use in source and binary forms, with or without
|
|
|
|
|
+ * modification, are permitted provided that the following conditions
|
|
|
|
|
+ * are met:
|
|
|
|
|
+ *
|
|
|
|
|
+ * 1. Redistributions of source code must retain the above copyright
|
|
|
|
|
+ * notice, this list of conditions and the following disclaimer.
|
|
|
|
|
+ * 2. Redistributions in binary form must reproduce the above
|
|
|
|
|
+ * copyright notice, this list of conditions and the following
|
|
|
|
|
+ * disclaimer in the documentation and/or other materials
|
|
|
|
|
+ * provided with the distribution.
|
|
|
|
|
+ *
|
|
|
|
|
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND
|
|
|
|
|
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
|
|
|
|
|
+ * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
|
|
|
|
|
+ * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS
|
|
|
|
|
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
|
|
|
|
|
+ * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
|
|
|
|
|
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
|
|
|
|
|
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
|
|
|
|
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
|
|
|
|
|
+ * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
|
|
|
|
|
+ * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
|
|
|
|
|
+ * OF SUCH DAMAGE.
|
|
|
|
|
+ *
|
|
|
|
|
+ * Reimplement by Jannson (20161018): compatible for golang version of https://github.com/klauspost/reedsolomon
|
|
|
|
|
+ */
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * The following parameter defines how many bits are used for
|
|
|
|
|
+ * field elements. The code supports any value from 2 to 16
|
|
|
|
|
+ * but fastest operation is achieved with 8 bit elements
|
|
|
|
|
+ * This is the only parameter you may want to change.
|
|
|
|
|
+ */
|
|
|
|
|
+#define GF_BITS 8 /* code over GF(2**GF_BITS) - change to suit */
|
|
|
|
|
+
|
|
|
|
|
+#include <stdio.h>
|
|
|
|
|
+#include <stdlib.h>
|
|
|
|
|
+#include <string.h>
|
|
|
|
|
+
|
|
|
|
|
+#include <assert.h>
|
|
|
|
|
+#include "rs.h"
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * stuff used for testing purposes only
|
|
|
|
|
+ */
|
|
|
|
|
+
|
|
|
|
|
+#ifdef TEST
|
|
|
|
|
+#define DEB(x)
|
|
|
|
|
+#define DDB(x) x
|
|
|
|
|
+#define DEBUG 0 /* minimal debugging */
|
|
|
|
|
+
|
|
|
|
|
+#include <sys/time.h>
|
|
|
|
|
+#define DIFF_T(a,b) \
|
|
|
|
|
+ (1+ 1000000*(a.tv_sec - b.tv_sec) + (a.tv_usec - b.tv_usec) )
|
|
|
|
|
+
|
|
|
|
|
+#define TICK(t) \
|
|
|
|
|
+ {struct timeval x ; \
|
|
|
|
|
+ gettimeofday(&x, NULL) ; \
|
|
|
|
|
+ t = x.tv_usec + 1000000* (x.tv_sec & 0xff ) ; \
|
|
|
|
|
+ }
|
|
|
|
|
+#define TOCK(t) \
|
|
|
|
|
+ { u_long t1 ; TICK(t1) ; \
|
|
|
|
|
+ if (t1 < t) t = 256000000 + t1 - t ; \
|
|
|
|
|
+ else t = t1 - t ; \
|
|
|
|
|
+ if (t == 0) t = 1 ;}
|
|
|
|
|
+
|
|
|
|
|
+u_long ticks[10]; /* vars for timekeeping */
|
|
|
|
|
+#else
|
|
|
|
|
+#define DEB(x)
|
|
|
|
|
+#define DDB(x)
|
|
|
|
|
+#define TICK(x)
|
|
|
|
|
+#define TOCK(x)
|
|
|
|
|
+#endif /* TEST */
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * You should not need to change anything beyond this point.
|
|
|
|
|
+ * The first part of the file implements linear algebra in GF.
|
|
|
|
|
+ *
|
|
|
|
|
+ * gf is the type used to store an element of the Galois Field.
|
|
|
|
|
+ * Must constain at least GF_BITS bits.
|
|
|
|
|
+ *
|
|
|
|
|
+ * Note: unsigned char will work up to GF(256) but int seems to run
|
|
|
|
|
+ * faster on the Pentium. We use int whenever have to deal with an
|
|
|
|
|
+ * index, since they are generally faster.
|
|
|
|
|
+ */
|
|
|
|
|
+/*
|
|
|
|
|
+ * AK: Udpcast only uses GF_BITS=8. Remove other possibilities
|
|
|
|
|
+ */
|
|
|
|
|
+#if (GF_BITS != 8)
|
|
|
|
|
+#error "GF_BITS must be 8"
|
|
|
|
|
+#endif
|
|
|
|
|
+typedef unsigned char gf;
|
|
|
|
|
+
|
|
|
|
|
+#define GF_SIZE ((1 << GF_BITS) - 1) /* powers of \alpha */
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * Primitive polynomials - see Lin & Costello, Appendix A,
|
|
|
|
|
+ * and Lee & Messerschmitt, p. 453.
|
|
|
|
|
+ */
|
|
|
|
|
+static char *allPp[] = { /* GF_BITS polynomial */
|
|
|
|
|
+ NULL, /* 0 no code */
|
|
|
|
|
+ NULL, /* 1 no code */
|
|
|
|
|
+ "111", /* 2 1+x+x^2 */
|
|
|
|
|
+ "1101", /* 3 1+x+x^3 */
|
|
|
|
|
+ "11001", /* 4 1+x+x^4 */
|
|
|
|
|
+ "101001", /* 5 1+x^2+x^5 */
|
|
|
|
|
+ "1100001", /* 6 1+x+x^6 */
|
|
|
|
|
+ "10010001", /* 7 1 + x^3 + x^7 */
|
|
|
|
|
+ "101110001", /* 8 1+x^2+x^3+x^4+x^8 */
|
|
|
|
|
+ "1000100001", /* 9 1+x^4+x^9 */
|
|
|
|
|
+ "10010000001", /* 10 1+x^3+x^10 */
|
|
|
|
|
+ "101000000001", /* 11 1+x^2+x^11 */
|
|
|
|
|
+ "1100101000001", /* 12 1+x+x^4+x^6+x^12 */
|
|
|
|
|
+ "11011000000001", /* 13 1+x+x^3+x^4+x^13 */
|
|
|
|
|
+ "110000100010001", /* 14 1+x+x^6+x^10+x^14 */
|
|
|
|
|
+ "1100000000000001", /* 15 1+x+x^15 */
|
|
|
|
|
+ "11010000000010001" /* 16 1+x+x^3+x^12+x^16 */
|
|
|
|
|
+};
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * To speed up computations, we have tables for logarithm, exponent
|
|
|
|
|
+ * and inverse of a number. If GF_BITS <= 8, we use a table for
|
|
|
|
|
+ * multiplication as well (it takes 64K, no big deal even on a PDA,
|
|
|
|
|
+ * especially because it can be pre-initialized an put into a ROM!),
|
|
|
|
|
+ * otherwhise we use a table of logarithms.
|
|
|
|
|
+ * In any case the macro gf_mul(x,y) takes care of multiplications.
|
|
|
|
|
+ */
|
|
|
|
|
+
|
|
|
|
|
+static gf gf_exp[2*GF_SIZE]; /* index->poly form conversion table */
|
|
|
|
|
+static int gf_log[GF_SIZE + 1]; /* Poly->index form conversion table */
|
|
|
|
|
+static gf inverse[GF_SIZE+1]; /* inverse of field elem. */
|
|
|
|
|
+/* inv[\alpha**i]=\alpha**(GF_SIZE-i-1) */
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * modnn(x) computes x % GF_SIZE, where GF_SIZE is 2**GF_BITS - 1,
|
|
|
|
|
+ * without a slow divide.
|
|
|
|
|
+ */
|
|
|
|
|
+static inline gf
|
|
|
|
|
+modnn(int x)
|
|
|
|
|
+{
|
|
|
|
|
+ while (x >= GF_SIZE) {
|
|
|
|
|
+ x -= GF_SIZE;
|
|
|
|
|
+ x = (x >> GF_BITS) + (x & GF_SIZE);
|
|
|
|
|
+ }
|
|
|
|
|
+ return x;
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+#define SWAP(a,b,t) {t tmp; tmp=a; a=b; b=tmp;}
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * gf_mul(x,y) multiplies two numbers. If GF_BITS<=8, it is much
|
|
|
|
|
+ * faster to use a multiplication table.
|
|
|
|
|
+ *
|
|
|
|
|
+ * USE_GF_MULC, GF_MULC0(c) and GF_ADDMULC(x) can be used when multiplying
|
|
|
|
|
+ * many numbers by the same constant. In this case the first
|
|
|
|
|
+ * call sets the constant, and others perform the multiplications.
|
|
|
|
|
+ * A value related to the multiplication is held in a local variable
|
|
|
|
|
+ * declared with USE_GF_MULC . See usage in addmul1().
|
|
|
|
|
+ */
|
|
|
|
|
+#ifdef _MSC_VER
|
|
|
|
|
+__declspec(align(16))
|
|
|
|
|
+#else
|
|
|
|
|
+_Alignas(16)
|
|
|
|
|
+#endif
|
|
|
|
|
+static gf gf_mul_table[(GF_SIZE + 1) * (GF_SIZE + 1)];
|
|
|
|
|
+
|
|
|
|
|
+#define gf_mul(x,y) gf_mul_table[(x<<8)+y]
|
|
|
|
|
+
|
|
|
|
|
+#define USE_GF_MULC register gf * __gf_mulc_
|
|
|
|
|
+#define GF_MULC0(c) __gf_mulc_ = &gf_mul_table[(c)<<8]
|
|
|
|
|
+#define GF_ADDMULC(dst, x) dst ^= __gf_mulc_[x]
|
|
|
|
|
+#define GF_MULC(dst, x) dst = __gf_mulc_[x]
|
|
|
|
|
+
|
|
|
|
|
+static void
|
|
|
|
|
+init_mul_table(void)
|
|
|
|
|
+{
|
|
|
|
|
+ int i, j;
|
|
|
|
|
+ for (i=0; i< GF_SIZE+1; i++)
|
|
|
|
|
+ for (j=0; j< GF_SIZE+1; j++)
|
|
|
|
|
+ gf_mul_table[(i<<8)+j] = gf_exp[modnn(gf_log[i] + gf_log[j]) ] ;
|
|
|
|
|
+
|
|
|
|
|
+ for (j=0; j< GF_SIZE+1; j++)
|
|
|
|
|
+ gf_mul_table[j] = gf_mul_table[j<<8] = 0;
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * Generate GF(2**m) from the irreducible polynomial p(X) in p[0]..p[m]
|
|
|
|
|
+ * Lookup tables:
|
|
|
|
|
+ * index->polynomial form gf_exp[] contains j= \alpha^i;
|
|
|
|
|
+ * polynomial form -> index form gf_log[ j = \alpha^i ] = i
|
|
|
|
|
+ * \alpha=x is the primitive element of GF(2^m)
|
|
|
|
|
+ *
|
|
|
|
|
+ * For efficiency, gf_exp[] has size 2*GF_SIZE, so that a simple
|
|
|
|
|
+ * multiplication of two numbers can be resolved without calling modnn
|
|
|
|
|
+ */
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * initialize the data structures used for computations in GF.
|
|
|
|
|
+ */
|
|
|
|
|
+static void
|
|
|
|
|
+generate_gf(void)
|
|
|
|
|
+{
|
|
|
|
|
+ int i;
|
|
|
|
|
+ gf mask;
|
|
|
|
|
+ char *Pp = allPp[GF_BITS] ;
|
|
|
|
|
+
|
|
|
|
|
+ mask = 1; /* x ** 0 = 1 */
|
|
|
|
|
+ gf_exp[GF_BITS] = 0; /* will be updated at the end of the 1st loop */
|
|
|
|
|
+ /*
|
|
|
|
|
+ * first, generate the (polynomial representation of) powers of \alpha,
|
|
|
|
|
+ * which are stored in gf_exp[i] = \alpha ** i .
|
|
|
|
|
+ * At the same time build gf_log[gf_exp[i]] = i .
|
|
|
|
|
+ * The first GF_BITS powers are simply bits shifted to the left.
|
|
|
|
|
+ */
|
|
|
|
|
+ for (i = 0; i < GF_BITS; i++, mask <<= 1 ) {
|
|
|
|
|
+ gf_exp[i] = mask;
|
|
|
|
|
+ gf_log[gf_exp[i]] = i;
|
|
|
|
|
+ /*
|
|
|
|
|
+ * If Pp[i] == 1 then \alpha ** i occurs in poly-repr
|
|
|
|
|
+ * gf_exp[GF_BITS] = \alpha ** GF_BITS
|
|
|
|
|
+ */
|
|
|
|
|
+ if ( Pp[i] == '1' )
|
|
|
|
|
+ gf_exp[GF_BITS] ^= mask;
|
|
|
|
|
+ }
|
|
|
|
|
+ /*
|
|
|
|
|
+ * now gf_exp[GF_BITS] = \alpha ** GF_BITS is complete, so can als
|
|
|
|
|
+ * compute its inverse.
|
|
|
|
|
+ */
|
|
|
|
|
+ gf_log[gf_exp[GF_BITS]] = GF_BITS;
|
|
|
|
|
+ /*
|
|
|
|
|
+ * Poly-repr of \alpha ** (i+1) is given by poly-repr of
|
|
|
|
|
+ * \alpha ** i shifted left one-bit and accounting for any
|
|
|
|
|
+ * \alpha ** GF_BITS term that may occur when poly-repr of
|
|
|
|
|
+ * \alpha ** i is shifted.
|
|
|
|
|
+ */
|
|
|
|
|
+ mask = 1 << (GF_BITS - 1 ) ;
|
|
|
|
|
+ for (i = GF_BITS + 1; i < GF_SIZE; i++) {
|
|
|
|
|
+ if (gf_exp[i - 1] >= mask)
|
|
|
|
|
+ gf_exp[i] = gf_exp[GF_BITS] ^ ((gf_exp[i - 1] ^ mask) << 1);
|
|
|
|
|
+ else
|
|
|
|
|
+ gf_exp[i] = gf_exp[i - 1] << 1;
|
|
|
|
|
+ gf_log[gf_exp[i]] = i;
|
|
|
|
|
+ }
|
|
|
|
|
+ /*
|
|
|
|
|
+ * log(0) is not defined, so use a special value
|
|
|
|
|
+ */
|
|
|
|
|
+ gf_log[0] = GF_SIZE ;
|
|
|
|
|
+ /* set the extended gf_exp values for fast multiply */
|
|
|
|
|
+ for (i = 0 ; i < GF_SIZE ; i++)
|
|
|
|
|
+ gf_exp[i + GF_SIZE] = gf_exp[i] ;
|
|
|
|
|
+
|
|
|
|
|
+ /*
|
|
|
|
|
+ * again special cases. 0 has no inverse. This used to
|
|
|
|
|
+ * be initialized to GF_SIZE, but it should make no difference
|
|
|
|
|
+ * since noone is supposed to read from here.
|
|
|
|
|
+ */
|
|
|
|
|
+ inverse[0] = 0 ;
|
|
|
|
|
+ inverse[1] = 1;
|
|
|
|
|
+ for (i=2; i<=GF_SIZE; i++)
|
|
|
|
|
+ inverse[i] = gf_exp[GF_SIZE-gf_log[i]];
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * Various linear algebra operations that i use often.
|
|
|
|
|
+ */
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * addmul() computes dst[] = dst[] + c * src[]
|
|
|
|
|
+ * This is used often, so better optimize it! Currently the loop is
|
|
|
|
|
+ * unrolled 16 times, a good value for 486 and pentium-class machines.
|
|
|
|
|
+ * The case c=0 is also optimized, whereas c=1 is not. These
|
|
|
|
|
+ * calls are unfrequent in my typical apps so I did not bother.
|
|
|
|
|
+ *
|
|
|
|
|
+ * Note that gcc on
|
|
|
|
|
+ */
|
|
|
|
|
+#if 0
|
|
|
|
|
+#define addmul(dst, src, c, sz) \
|
|
|
|
|
+ if (c != 0) addmul1(dst, src, c, sz)
|
|
|
|
|
+#endif
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+#define UNROLL 16 /* 1, 4, 8, 16 */
|
|
|
|
|
+static void
|
|
|
|
|
+slow_addmul1(gf *dst1, gf *src1, gf c, int sz)
|
|
|
|
|
+{
|
|
|
|
|
+ USE_GF_MULC ;
|
|
|
|
|
+ register gf *dst = dst1, *src = src1 ;
|
|
|
|
|
+ gf *lim = &dst[sz - UNROLL + 1] ;
|
|
|
|
|
+
|
|
|
|
|
+ GF_MULC0(c) ;
|
|
|
|
|
+
|
|
|
|
|
+#if (UNROLL > 1) /* unrolling by 8/16 is quite effective on the pentium */
|
|
|
|
|
+ for (; dst < lim ; dst += UNROLL, src += UNROLL ) {
|
|
|
|
|
+ GF_ADDMULC( dst[0] , src[0] );
|
|
|
|
|
+ GF_ADDMULC( dst[1] , src[1] );
|
|
|
|
|
+ GF_ADDMULC( dst[2] , src[2] );
|
|
|
|
|
+ GF_ADDMULC( dst[3] , src[3] );
|
|
|
|
|
+#if (UNROLL > 4)
|
|
|
|
|
+ GF_ADDMULC( dst[4] , src[4] );
|
|
|
|
|
+ GF_ADDMULC( dst[5] , src[5] );
|
|
|
|
|
+ GF_ADDMULC( dst[6] , src[6] );
|
|
|
|
|
+ GF_ADDMULC( dst[7] , src[7] );
|
|
|
|
|
+#endif
|
|
|
|
|
+#if (UNROLL > 8)
|
|
|
|
|
+ GF_ADDMULC( dst[8] , src[8] );
|
|
|
|
|
+ GF_ADDMULC( dst[9] , src[9] );
|
|
|
|
|
+ GF_ADDMULC( dst[10] , src[10] );
|
|
|
|
|
+ GF_ADDMULC( dst[11] , src[11] );
|
|
|
|
|
+ GF_ADDMULC( dst[12] , src[12] );
|
|
|
|
|
+ GF_ADDMULC( dst[13] , src[13] );
|
|
|
|
|
+ GF_ADDMULC( dst[14] , src[14] );
|
|
|
|
|
+ GF_ADDMULC( dst[15] , src[15] );
|
|
|
|
|
+#endif
|
|
|
|
|
+ }
|
|
|
|
|
+#endif
|
|
|
|
|
+ lim += UNROLL - 1 ;
|
|
|
|
|
+ for (; dst < lim; dst++, src++ ) /* final components */
|
|
|
|
|
+ GF_ADDMULC( *dst , *src );
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+# define addmul1 slow_addmul1
|
|
|
|
|
+
|
|
|
|
|
+static void addmul(gf *dst, gf *src, gf c, int sz) {
|
|
|
|
|
+ // fprintf(stderr, "Dst=%p Src=%p, gf=%02x sz=%d\n", dst, src, c, sz);
|
|
|
|
|
+ if (c != 0) addmul1(dst, src, c, sz);
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * mul() computes dst[] = c * src[]
|
|
|
|
|
+ * This is used often, so better optimize it! Currently the loop is
|
|
|
|
|
+ * unrolled 16 times, a good value for 486 and pentium-class machines.
|
|
|
|
|
+ * The case c=0 is also optimized, whereas c=1 is not. These
|
|
|
|
|
+ * calls are unfrequent in my typical apps so I did not bother.
|
|
|
|
|
+ *
|
|
|
|
|
+ * Note that gcc on
|
|
|
|
|
+ */
|
|
|
|
|
+#if 0
|
|
|
|
|
+#define mul(dst, src, c, sz) \
|
|
|
|
|
+ do { if (c != 0) mul1(dst, src, c, sz); else memset(dst, 0, c); } while(0)
|
|
|
|
|
+#endif
|
|
|
|
|
+
|
|
|
|
|
+#define UNROLL 16 /* 1, 4, 8, 16 */
|
|
|
|
|
+static void
|
|
|
|
|
+slow_mul1(gf *dst1, gf *src1, gf c, int sz)
|
|
|
|
|
+{
|
|
|
|
|
+ USE_GF_MULC ;
|
|
|
|
|
+ register gf *dst = dst1, *src = src1 ;
|
|
|
|
|
+ gf *lim = &dst[sz - UNROLL + 1] ;
|
|
|
|
|
+
|
|
|
|
|
+ GF_MULC0(c) ;
|
|
|
|
|
+
|
|
|
|
|
+#if (UNROLL > 1) /* unrolling by 8/16 is quite effective on the pentium */
|
|
|
|
|
+ for (; dst < lim ; dst += UNROLL, src += UNROLL ) {
|
|
|
|
|
+ GF_MULC( dst[0] , src[0] );
|
|
|
|
|
+ GF_MULC( dst[1] , src[1] );
|
|
|
|
|
+ GF_MULC( dst[2] , src[2] );
|
|
|
|
|
+ GF_MULC( dst[3] , src[3] );
|
|
|
|
|
+#if (UNROLL > 4)
|
|
|
|
|
+ GF_MULC( dst[4] , src[4] );
|
|
|
|
|
+ GF_MULC( dst[5] , src[5] );
|
|
|
|
|
+ GF_MULC( dst[6] , src[6] );
|
|
|
|
|
+ GF_MULC( dst[7] , src[7] );
|
|
|
|
|
+#endif
|
|
|
|
|
+#if (UNROLL > 8)
|
|
|
|
|
+ GF_MULC( dst[8] , src[8] );
|
|
|
|
|
+ GF_MULC( dst[9] , src[9] );
|
|
|
|
|
+ GF_MULC( dst[10] , src[10] );
|
|
|
|
|
+ GF_MULC( dst[11] , src[11] );
|
|
|
|
|
+ GF_MULC( dst[12] , src[12] );
|
|
|
|
|
+ GF_MULC( dst[13] , src[13] );
|
|
|
|
|
+ GF_MULC( dst[14] , src[14] );
|
|
|
|
|
+ GF_MULC( dst[15] , src[15] );
|
|
|
|
|
+#endif
|
|
|
|
|
+ }
|
|
|
|
|
+#endif
|
|
|
|
|
+ lim += UNROLL - 1 ;
|
|
|
|
|
+ for (; dst < lim; dst++, src++ ) /* final components */
|
|
|
|
|
+ GF_MULC( *dst , *src );
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+# define mul1 slow_mul1
|
|
|
|
|
+
|
|
|
|
|
+static inline void mul(gf *dst, gf *src, gf c, int sz) {
|
|
|
|
|
+ /*fprintf(stderr, "%p = %02x * %p\n", dst, c, src);*/
|
|
|
|
|
+ if (c != 0) mul1(dst, src, c, sz); else memset(dst, 0, c);
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * invert_mat() takes a matrix and produces its inverse
|
|
|
|
|
+ * k is the size of the matrix.
|
|
|
|
|
+ * (Gauss-Jordan, adapted from Numerical Recipes in C)
|
|
|
|
|
+ * Return non-zero if singular.
|
|
|
|
|
+ */
|
|
|
|
|
+DEB( int pivloops=0; int pivswaps=0 ; /* diagnostic */)
|
|
|
|
|
+static int
|
|
|
|
|
+invert_mat(gf *src, int k)
|
|
|
|
|
+{
|
|
|
|
|
+ gf c, *p ;
|
|
|
|
|
+ int irow, icol, row, col, i, ix ;
|
|
|
|
|
+
|
|
|
|
|
+ int error = 1 ;
|
|
|
|
|
+ int *indxc = malloc(k*sizeof(int));
|
|
|
|
|
+ int *indxr = malloc(k*sizeof(int));
|
|
|
|
|
+ int *ipiv = malloc(k*sizeof(int));
|
|
|
|
|
+ gf *id_row = malloc(k*sizeof(gf));
|
|
|
|
|
+// int indxc[k];
|
|
|
|
|
+// int indxr[k];
|
|
|
|
|
+// int ipiv[k];
|
|
|
|
|
+// gf id_row[k];
|
|
|
|
|
+
|
|
|
|
|
+ memset(id_row, 0, k*sizeof(gf));
|
|
|
|
|
+ DEB( pivloops=0; pivswaps=0 ; /* diagnostic */ )
|
|
|
|
|
+ /*
|
|
|
|
|
+ * ipiv marks elements already used as pivots.
|
|
|
|
|
+ */
|
|
|
|
|
+ for (i = 0; i < k ; i++)
|
|
|
|
|
+ ipiv[i] = 0 ;
|
|
|
|
|
+
|
|
|
|
|
+ for (col = 0; col < k ; col++) {
|
|
|
|
|
+ gf *pivot_row ;
|
|
|
|
|
+ /*
|
|
|
|
|
+ * Zeroing column 'col', look for a non-zero element.
|
|
|
|
|
+ * First try on the diagonal, if it fails, look elsewhere.
|
|
|
|
|
+ */
|
|
|
|
|
+ irow = icol = -1 ;
|
|
|
|
|
+ if (ipiv[col] != 1 && src[col*k + col] != 0) {
|
|
|
|
|
+ irow = col ;
|
|
|
|
|
+ icol = col ;
|
|
|
|
|
+ goto found_piv ;
|
|
|
|
|
+ }
|
|
|
|
|
+ for (row = 0 ; row < k ; row++) {
|
|
|
|
|
+ if (ipiv[row] != 1) {
|
|
|
|
|
+ for (ix = 0 ; ix < k ; ix++) {
|
|
|
|
|
+ DEB( pivloops++ ; )
|
|
|
|
|
+ if (ipiv[ix] == 0) {
|
|
|
|
|
+ if (src[row*k + ix] != 0) {
|
|
|
|
|
+ irow = row ;
|
|
|
|
|
+ icol = ix ;
|
|
|
|
|
+ goto found_piv ;
|
|
|
|
|
+ }
|
|
|
|
|
+ } else if (ipiv[ix] > 1) {
|
|
|
|
|
+ fprintf(stderr, "singular matrix\n");
|
|
|
|
|
+ goto fail ;
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ if (icol == -1) {
|
|
|
|
|
+ fprintf(stderr, "XXX pivot not found!\n");
|
|
|
|
|
+ goto fail ;
|
|
|
|
|
+ }
|
|
|
|
|
+ found_piv:
|
|
|
|
|
+ ++(ipiv[icol]) ;
|
|
|
|
|
+ /*
|
|
|
|
|
+ * swap rows irow and icol, so afterwards the diagonal
|
|
|
|
|
+ * element will be correct. Rarely done, not worth
|
|
|
|
|
+ * optimizing.
|
|
|
|
|
+ */
|
|
|
|
|
+ if (irow != icol) {
|
|
|
|
|
+ for (ix = 0 ; ix < k ; ix++ ) {
|
|
|
|
|
+ SWAP( src[irow*k + ix], src[icol*k + ix], gf) ;
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ indxr[col] = irow ;
|
|
|
|
|
+ indxc[col] = icol ;
|
|
|
|
|
+ pivot_row = &src[icol*k] ;
|
|
|
|
|
+ c = pivot_row[icol] ;
|
|
|
|
|
+ if (c == 0) {
|
|
|
|
|
+ fprintf(stderr, "singular matrix 2\n");
|
|
|
|
|
+ goto fail ;
|
|
|
|
|
+ }
|
|
|
|
|
+ if (c != 1 ) { /* otherwhise this is a NOP */
|
|
|
|
|
+ /*
|
|
|
|
|
+ * this is done often , but optimizing is not so
|
|
|
|
|
+ * fruitful, at least in the obvious ways (unrolling)
|
|
|
|
|
+ */
|
|
|
|
|
+ DEB( pivswaps++ ; )
|
|
|
|
|
+ c = inverse[ c ] ;
|
|
|
|
|
+ pivot_row[icol] = 1 ;
|
|
|
|
|
+ for (ix = 0 ; ix < k ; ix++ )
|
|
|
|
|
+ pivot_row[ix] = gf_mul(c, pivot_row[ix] );
|
|
|
|
|
+ }
|
|
|
|
|
+ /*
|
|
|
|
|
+ * from all rows, remove multiples of the selected row
|
|
|
|
|
+ * to zero the relevant entry (in fact, the entry is not zero
|
|
|
|
|
+ * because we know it must be zero).
|
|
|
|
|
+ * (Here, if we know that the pivot_row is the identity,
|
|
|
|
|
+ * we can optimize the addmul).
|
|
|
|
|
+ */
|
|
|
|
|
+ id_row[icol] = 1;
|
|
|
|
|
+ if (memcmp(pivot_row, id_row, k*sizeof(gf)) != 0) {
|
|
|
|
|
+ for (p = src, ix = 0 ; ix < k ; ix++, p += k ) {
|
|
|
|
|
+ if (ix != icol) {
|
|
|
|
|
+ c = p[icol] ;
|
|
|
|
|
+ p[icol] = 0 ;
|
|
|
|
|
+ addmul(p, pivot_row, c, k );
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ id_row[icol] = 0;
|
|
|
|
|
+ } /* done all columns */
|
|
|
|
|
+ for (col = k-1 ; col >= 0 ; col-- ) {
|
|
|
|
|
+ if (indxr[col] <0 || indxr[col] >= k)
|
|
|
|
|
+ fprintf(stderr, "AARGH, indxr[col] %d\n", indxr[col]);
|
|
|
|
|
+ else if (indxc[col] <0 || indxc[col] >= k)
|
|
|
|
|
+ fprintf(stderr, "AARGH, indxc[col] %d\n", indxc[col]);
|
|
|
|
|
+ else
|
|
|
|
|
+ if (indxr[col] != indxc[col] ) {
|
|
|
|
|
+ for (row = 0 ; row < k ; row++ ) {
|
|
|
|
|
+ SWAP( src[row*k + indxr[col]], src[row*k + indxc[col]], gf) ;
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ error = 0 ;
|
|
|
|
|
+ fail:
|
|
|
|
|
+ free(indxc);
|
|
|
|
|
+ free(indxr);
|
|
|
|
|
+ free(ipiv);
|
|
|
|
|
+ free(id_row);
|
|
|
|
|
+ return error ;
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+static int fec_initialized = 0 ;
|
|
|
|
|
+
|
|
|
|
|
+void fec_init(void)
|
|
|
|
|
+{
|
|
|
|
|
+ TICK(ticks[0]);
|
|
|
|
|
+ generate_gf();
|
|
|
|
|
+ TOCK(ticks[0]);
|
|
|
|
|
+ DDB(fprintf(stderr, "generate_gf took %ldus\n", ticks[0]);)
|
|
|
|
|
+ TICK(ticks[0]);
|
|
|
|
|
+ init_mul_table();
|
|
|
|
|
+ TOCK(ticks[0]);
|
|
|
|
|
+ DDB(fprintf(stderr, "init_mul_table took %ldus\n", ticks[0]);)
|
|
|
|
|
+ fec_initialized = 1 ;
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+#ifdef PROFILE
|
|
|
|
|
+#ifdef __x86_64__
|
|
|
|
|
+static long long rdtsc(void)
|
|
|
|
|
+{
|
|
|
|
|
+ unsigned long low, hi;
|
|
|
|
|
+ asm volatile ("rdtsc" : "=d" (hi), "=a" (low));
|
|
|
|
|
+ return ( (((long long)hi) << 32) | ((long long) low));
|
|
|
|
|
+}
|
|
|
|
|
+#elif __arm__
|
|
|
|
|
+static long long rdtsc(void)
|
|
|
|
|
+{
|
|
|
|
|
+ u64 val;
|
|
|
|
|
+ asm volatile("mrs %0, cntvct_el0" : "=r" (val));
|
|
|
|
|
+ return val;
|
|
|
|
|
+}
|
|
|
|
|
+#endif
|
|
|
|
|
+
|
|
|
|
|
+void print_matrix1(gf* matrix, int nrows, int ncols) {
|
|
|
|
|
+ int i, j;
|
|
|
|
|
+ printf("matrix (%d,%d):\n", nrows, ncols);
|
|
|
|
|
+ for(i = 0; i < nrows; i++) {
|
|
|
|
|
+ for(j = 0; j < ncols; j++) {
|
|
|
|
|
+ printf("%6d ", matrix[i*ncols + j]);
|
|
|
|
|
+ }
|
|
|
|
|
+ printf("\n");
|
|
|
|
|
+ }
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+void print_matrix2(gf** matrix, int nrows, int ncols) {
|
|
|
|
|
+ int i, j;
|
|
|
|
|
+ printf("matrix (%d,%d):\n", nrows, ncols);
|
|
|
|
|
+ for(i = 0; i < nrows; i++) {
|
|
|
|
|
+ for(j = 0; j < ncols; j++) {
|
|
|
|
|
+ printf("%6d ", matrix[i][j]);
|
|
|
|
|
+ }
|
|
|
|
|
+ printf("\n");
|
|
|
|
|
+ }
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+#endif
|
|
|
|
|
+
|
|
|
|
|
+/* y = a**n */
|
|
|
|
|
+static gf galExp(gf a, gf n) {
|
|
|
|
|
+ int logA;
|
|
|
|
|
+ int logResult;
|
|
|
|
|
+ if(0 == n) {
|
|
|
|
|
+ return 1;
|
|
|
|
|
+ }
|
|
|
|
|
+ if(0 == a) {
|
|
|
|
|
+ return 0;
|
|
|
|
|
+ }
|
|
|
|
|
+ logA = gf_log[a];
|
|
|
|
|
+ logResult = logA * n;
|
|
|
|
|
+ while(logResult >= 255) {
|
|
|
|
|
+ logResult -= 255;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ return gf_exp[logResult];
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+static inline gf galMultiply(gf a, gf b) {
|
|
|
|
|
+ return gf_mul_table[ ((int)a << 8) + (int)b ];
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+static gf* vandermonde(int nrows, int ncols) {
|
|
|
|
|
+ int row, col, ptr;
|
|
|
|
|
+ gf* matrix = (gf*)RS_MALLOC(nrows * ncols);
|
|
|
|
|
+ if(NULL != matrix) {
|
|
|
|
|
+ ptr = 0;
|
|
|
|
|
+ for(row = 0; row < nrows; row++) {
|
|
|
|
|
+ for(col = 0; col < ncols; col++) {
|
|
|
|
|
+ matrix[ptr++] = galExp((gf)row, (gf)col);
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ return matrix;
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/*
|
|
|
|
|
+ * Not check for input params
|
|
|
|
|
+ * */
|
|
|
|
|
+static gf* sub_matrix(gf* matrix, int rmin, int cmin, int rmax, int cmax, int nrows, int ncols) {
|
|
|
|
|
+ int i, j, ptr = 0;
|
|
|
|
|
+ gf* new_m = (gf*)RS_MALLOC( (rmax-rmin) * (cmax-cmin) );
|
|
|
|
|
+ if(NULL != new_m) {
|
|
|
|
|
+ for(i = rmin; i < rmax; i++) {
|
|
|
|
|
+ for(j = cmin; j < cmax; j++) {
|
|
|
|
|
+ new_m[ptr++] = matrix[i*ncols + j];
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ return new_m;
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/* y = a.dot(b) */
|
|
|
|
|
+static gf* multiply1(gf *a, int ar, int ac, gf *b, int br, int bc) {
|
|
|
|
|
+ gf *new_m, tg;
|
|
|
|
|
+ int r, c, i, ptr = 0;
|
|
|
|
|
+
|
|
|
|
|
+ assert(ac == br);
|
|
|
|
|
+ new_m = (gf*)RS_CALLOC(1, ar*bc);
|
|
|
|
|
+ if(NULL != new_m) {
|
|
|
|
|
+
|
|
|
|
|
+ /* this multiply is slow */
|
|
|
|
|
+ for(r = 0; r < ar; r++) {
|
|
|
|
|
+ for(c = 0; c < bc; c++) {
|
|
|
|
|
+ tg = 0;
|
|
|
|
|
+ for(i = 0; i < ac; i++) {
|
|
|
|
|
+ /* tg ^= gf_mul_table[ ((int)a[r*ac+i] << 8) + (int)b[i*bc+c] ]; */
|
|
|
|
|
+ tg ^= galMultiply(a[r*ac+i], b[i*bc+c]);
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ new_m[ptr++] = tg;
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ return new_m;
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/* copy from golang rs version */
|
|
|
|
|
+static inline int code_some_shards(gf* matrixRows, gf** inputs, gf** outputs,
|
|
|
|
|
+ int dataShards, int outputCount, int byteCount) {
|
|
|
|
|
+ gf* in;
|
|
|
|
|
+ int iRow, c;
|
|
|
|
|
+ for(c = 0; c < dataShards; c++) {
|
|
|
|
|
+ in = inputs[c];
|
|
|
|
|
+ for(iRow = 0; iRow < outputCount; iRow++) {
|
|
|
|
|
+ if(0 == c) {
|
|
|
|
|
+ mul(outputs[iRow], in, matrixRows[iRow*dataShards+c], byteCount);
|
|
|
|
|
+ } else {
|
|
|
|
|
+ addmul(outputs[iRow], in, matrixRows[iRow*dataShards+c], byteCount);
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ return 0;
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+reed_solomon* reed_solomon_new(int data_shards, int parity_shards) {
|
|
|
|
|
+ gf* vm = NULL;
|
|
|
|
|
+ gf* top = NULL;
|
|
|
|
|
+ int err = 0;
|
|
|
|
|
+ reed_solomon* rs = NULL;
|
|
|
|
|
+
|
|
|
|
|
+ /* MUST use fec_init once time first */
|
|
|
|
|
+ assert(fec_initialized);
|
|
|
|
|
+
|
|
|
|
|
+ do {
|
|
|
|
|
+ rs = (reed_solomon*) RS_MALLOC(sizeof(reed_solomon));
|
|
|
|
|
+ if(NULL == rs) {
|
|
|
|
|
+ return NULL;
|
|
|
|
|
+ }
|
|
|
|
|
+ rs->data_shards = data_shards;
|
|
|
|
|
+ rs->parity_shards = parity_shards;
|
|
|
|
|
+ rs->shards = (data_shards + parity_shards);
|
|
|
|
|
+ rs->m = NULL;
|
|
|
|
|
+ rs->parity = NULL;
|
|
|
|
|
+
|
|
|
|
|
+ if(rs->shards > DATA_SHARDS_MAX || data_shards <= 0 || parity_shards <= 0) {
|
|
|
|
|
+ err = 1;
|
|
|
|
|
+ break;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ vm = vandermonde(rs->shards, rs->data_shards);
|
|
|
|
|
+ if(NULL == vm) {
|
|
|
|
|
+ err = 2;
|
|
|
|
|
+ break;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ top = sub_matrix(vm, 0, 0, data_shards, data_shards, rs->shards, data_shards);
|
|
|
|
|
+ if(NULL == top) {
|
|
|
|
|
+ err = 3;
|
|
|
|
|
+ break;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ err = invert_mat(top, data_shards);
|
|
|
|
|
+ assert(0 == err);
|
|
|
|
|
+
|
|
|
|
|
+ rs->m = multiply1(vm, rs->shards, data_shards, top, data_shards, data_shards);
|
|
|
|
|
+ if(NULL == rs->m) {
|
|
|
|
|
+ err = 4;
|
|
|
|
|
+ break;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ rs->parity = sub_matrix(rs->m, data_shards, 0, rs->shards, data_shards, rs->shards, data_shards);
|
|
|
|
|
+ if(NULL == rs->parity) {
|
|
|
|
|
+ err = 5;
|
|
|
|
|
+ break;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ RS_FREE(vm);
|
|
|
|
|
+ RS_FREE(top);
|
|
|
|
|
+ vm = NULL;
|
|
|
|
|
+ top = NULL;
|
|
|
|
|
+ return rs;
|
|
|
|
|
+
|
|
|
|
|
+ } while(0);
|
|
|
|
|
+
|
|
|
|
|
+ fprintf(stderr, "err=%d\n", err);
|
|
|
|
|
+ if(NULL != vm) {
|
|
|
|
|
+ RS_FREE(vm);
|
|
|
|
|
+ }
|
|
|
|
|
+ if(NULL != top) {
|
|
|
|
|
+ RS_FREE(top);
|
|
|
|
|
+ }
|
|
|
|
|
+ if(NULL != rs) {
|
|
|
|
|
+ if(NULL != rs->m) {
|
|
|
|
|
+ RS_FREE(rs->m);
|
|
|
|
|
+ }
|
|
|
|
|
+ if(NULL != rs->parity) {
|
|
|
|
|
+ RS_FREE(rs->parity);
|
|
|
|
|
+ }
|
|
|
|
|
+ RS_FREE(rs);
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ return NULL;
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+void reed_solomon_release(reed_solomon* rs) {
|
|
|
|
|
+ if(NULL != rs) {
|
|
|
|
|
+ if(NULL != rs->m) {
|
|
|
|
|
+ RS_FREE(rs->m);
|
|
|
|
|
+ }
|
|
|
|
|
+ if(NULL != rs->parity) {
|
|
|
|
|
+ RS_FREE(rs->parity);
|
|
|
|
|
+ }
|
|
|
|
|
+ RS_FREE(rs);
|
|
|
|
|
+ }
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/**
|
|
|
|
|
+ * encode one shard
|
|
|
|
|
+ * input:
|
|
|
|
|
+ * rs
|
|
|
|
|
+ * data_blocks[rs->data_shards][block_size]
|
|
|
|
|
+ * fec_blocks[rs->data_shards][block_size]
|
|
|
|
|
+ * */
|
|
|
|
|
+int reed_solomon_encode(reed_solomon* rs,
|
|
|
|
|
+ unsigned char** data_blocks,
|
|
|
|
|
+ unsigned char** fec_blocks,
|
|
|
|
|
+ int block_size) {
|
|
|
|
|
+ assert(NULL != rs && NULL != rs->parity);
|
|
|
|
|
+
|
|
|
|
|
+ return code_some_shards(rs->parity, data_blocks, fec_blocks
|
|
|
|
|
+ , rs->data_shards, rs->parity_shards, block_size);
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/**
|
|
|
|
|
+ * decode one shard
|
|
|
|
|
+ * input:
|
|
|
|
|
+ * rs
|
|
|
|
|
+ * original data_blocks[rs->data_shards][block_size]
|
|
|
|
|
+ * dec_fec_blocks[nr_fec_blocks][block_size]
|
|
|
|
|
+ * fec_block_nos: fec pos number in original fec_blocks
|
|
|
|
|
+ * erased_blocks: erased blocks in original data_blocks
|
|
|
|
|
+ * nr_fec_blocks: the number of erased blocks
|
|
|
|
|
+ * */
|
|
|
|
|
+int reed_solomon_decode(reed_solomon* rs,
|
|
|
|
|
+ unsigned char **data_blocks,
|
|
|
|
|
+ int block_size,
|
|
|
|
|
+ unsigned char **dec_fec_blocks,
|
|
|
|
|
+ unsigned int *fec_block_nos,
|
|
|
|
|
+ unsigned int *erased_blocks,
|
|
|
|
|
+ int nr_fec_blocks) {
|
|
|
|
|
+ /* use stack instead of malloc, define a small number of DATA_SHARDS_MAX to save memory */
|
|
|
|
|
+ gf dataDecodeMatrix[DATA_SHARDS_MAX*DATA_SHARDS_MAX];
|
|
|
|
|
+ unsigned char* subShards[DATA_SHARDS_MAX];
|
|
|
|
|
+ unsigned char* outputs[DATA_SHARDS_MAX];
|
|
|
|
|
+ gf* m = rs->m;
|
|
|
|
|
+ int i, j, c, swap, subMatrixRow, dataShards, nos, nshards;
|
|
|
|
|
+
|
|
|
|
|
+ /* the erased_blocks should always sorted
|
|
|
|
|
+ * if sorted, nr_fec_blocks times to check it
|
|
|
|
|
+ * if not, sort it here
|
|
|
|
|
+ * */
|
|
|
|
|
+ for(i = 0; i < nr_fec_blocks; i++) {
|
|
|
|
|
+ swap = 0;
|
|
|
|
|
+ for(j = i+1; j < nr_fec_blocks; j++) {
|
|
|
|
|
+ if(erased_blocks[i] > erased_blocks[j]) {
|
|
|
|
|
+ /* the prefix is bigger than the following, swap */
|
|
|
|
|
+ c = erased_blocks[i];
|
|
|
|
|
+ erased_blocks[i] = erased_blocks[j];
|
|
|
|
|
+ erased_blocks[j] = c;
|
|
|
|
|
+
|
|
|
|
|
+ swap = 1;
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ //printf("swap:%d\n", swap);
|
|
|
|
|
+ if(!swap) {
|
|
|
|
|
+ //already sorted or sorted ok
|
|
|
|
|
+ break;
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ j = 0;
|
|
|
|
|
+ subMatrixRow = 0;
|
|
|
|
|
+ nos = 0;
|
|
|
|
|
+ nshards = 0;
|
|
|
|
|
+ dataShards = rs->data_shards;
|
|
|
|
|
+ for(i = 0; i < dataShards; i++) {
|
|
|
|
|
+ if(j < nr_fec_blocks && i == erased_blocks[j]) {
|
|
|
|
|
+ //ignore the invalid block
|
|
|
|
|
+ j++;
|
|
|
|
|
+ } else {
|
|
|
|
|
+ /* this row is ok */
|
|
|
|
|
+ for(c = 0; c < dataShards; c++) {
|
|
|
|
|
+ dataDecodeMatrix[subMatrixRow*dataShards + c] = m[i*dataShards + c];
|
|
|
|
|
+ }
|
|
|
|
|
+ subShards[subMatrixRow] = data_blocks[i];
|
|
|
|
|
+ subMatrixRow++;
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ for(i = 0; i < nr_fec_blocks && subMatrixRow < dataShards; i++) {
|
|
|
|
|
+ subShards[subMatrixRow] = dec_fec_blocks[i];
|
|
|
|
|
+ j = dataShards + fec_block_nos[i];
|
|
|
|
|
+ for(c = 0; c < dataShards; c++) {
|
|
|
|
|
+ dataDecodeMatrix[subMatrixRow*dataShards + c] = m[j*dataShards + c]; //use spefic pos of original fec_blocks
|
|
|
|
|
+ }
|
|
|
|
|
+ subMatrixRow++;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ if(subMatrixRow < dataShards) {
|
|
|
|
|
+ //cannot correct
|
|
|
|
|
+ return -1;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ invert_mat(dataDecodeMatrix, dataShards);
|
|
|
|
|
+ //printf("invert:\n");
|
|
|
|
|
+ //print_matrix1(dataDecodeMatrix, dataShards, dataShards);
|
|
|
|
|
+ //printf("nShards:\n");
|
|
|
|
|
+ //print_matrix2(subShards, dataShards, block_size);
|
|
|
|
|
+
|
|
|
|
|
+ for(i = 0; i < nr_fec_blocks; i++) {
|
|
|
|
|
+ j = erased_blocks[i];
|
|
|
|
|
+ outputs[i] = data_blocks[j];
|
|
|
|
|
+ //data_blocks[j][0] = 0;
|
|
|
|
|
+ memmove(dataDecodeMatrix+i*dataShards, dataDecodeMatrix+j*dataShards, dataShards);
|
|
|
|
|
+ }
|
|
|
|
|
+ //printf("subMatrixRow:\n");
|
|
|
|
|
+ //print_matrix1(dataDecodeMatrix, nr_fec_blocks, dataShards);
|
|
|
|
|
+
|
|
|
|
|
+ //printf("outputs:\n");
|
|
|
|
|
+ //print_matrix2(outputs, nr_fec_blocks, block_size);
|
|
|
|
|
+
|
|
|
|
|
+ return code_some_shards(dataDecodeMatrix, subShards, outputs,
|
|
|
|
|
+ dataShards, nr_fec_blocks, block_size);
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/**
|
|
|
|
|
+ * encode a big size of buffer
|
|
|
|
|
+ * input:
|
|
|
|
|
+ * rs
|
|
|
|
|
+ * nr_shards: assert(0 == nr_shards % rs->shards)
|
|
|
|
|
+ * shards[nr_shards][block_size]
|
|
|
|
|
+ * */
|
|
|
|
|
+int reed_solomon_encode2(reed_solomon* rs, unsigned char** shards, int nr_shards, int block_size) {
|
|
|
|
|
+ unsigned char** data_blocks;
|
|
|
|
|
+ unsigned char** fec_blocks;
|
|
|
|
|
+ int i, ds = rs->data_shards, ps = rs->parity_shards, ss = rs->shards;
|
|
|
|
|
+ i = nr_shards / ss;
|
|
|
|
|
+ data_blocks = shards;
|
|
|
|
|
+ fec_blocks = &shards[(i*ds)];
|
|
|
|
|
+
|
|
|
|
|
+ for(i = 0; i < nr_shards; i += ss) {
|
|
|
|
|
+ reed_solomon_encode(rs, data_blocks, fec_blocks, block_size);
|
|
|
|
|
+ data_blocks += ds;
|
|
|
|
|
+ fec_blocks += ps;
|
|
|
|
|
+ }
|
|
|
|
|
+ return 0;
|
|
|
|
|
+}
|
|
|
|
|
+
|
|
|
|
|
+/**
|
|
|
|
|
+ * reconstruct a big size of buffer
|
|
|
|
|
+ * input:
|
|
|
|
|
+ * rs
|
|
|
|
|
+ * nr_shards: assert(0 == nr_shards % rs->data_shards)
|
|
|
|
|
+ * shards[nr_shards][block_size]
|
|
|
|
|
+ * marks[nr_shards] marks as errors
|
|
|
|
|
+ * */
|
|
|
|
|
+int reed_solomon_reconstruct(reed_solomon* rs,
|
|
|
|
|
+ unsigned char** shards,
|
|
|
|
|
+ unsigned char* marks,
|
|
|
|
|
+ int nr_shards,
|
|
|
|
|
+ int block_size) {
|
|
|
|
|
+ unsigned char *dec_fec_blocks[DATA_SHARDS_MAX];
|
|
|
|
|
+ unsigned int fec_block_nos[DATA_SHARDS_MAX];
|
|
|
|
|
+ unsigned int erased_blocks[DATA_SHARDS_MAX];
|
|
|
|
|
+ unsigned char* fec_marks;
|
|
|
|
|
+ unsigned char **data_blocks, **fec_blocks;
|
|
|
|
|
+ int i, j, dn, pn, n;
|
|
|
|
|
+ int ds = rs->data_shards;
|
|
|
|
|
+ int ps = rs->parity_shards;
|
|
|
|
|
+ int err = 0;
|
|
|
|
|
+
|
|
|
|
|
+ data_blocks = shards;
|
|
|
|
|
+ n = nr_shards / rs->shards;
|
|
|
|
|
+ fec_marks = marks + n*ds; //after all data, is't fec marks
|
|
|
|
|
+ fec_blocks = shards + n*ds;
|
|
|
|
|
+
|
|
|
|
|
+ for(j = 0; j < n; j++) {
|
|
|
|
|
+ dn = 0;
|
|
|
|
|
+ for(i = 0; i < ds; i++) {
|
|
|
|
|
+ if(marks[i]) {
|
|
|
|
|
+ //errors
|
|
|
|
|
+ erased_blocks[dn++] = i;
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ if(dn > 0) {
|
|
|
|
|
+ pn = 0;
|
|
|
|
|
+ for(i = 0; i < ps && pn < dn; i++) {
|
|
|
|
|
+ if(!fec_marks[i]) {
|
|
|
|
|
+ //got valid fec row
|
|
|
|
|
+ fec_block_nos[pn] = i;
|
|
|
|
|
+ dec_fec_blocks[pn] = fec_blocks[i];
|
|
|
|
|
+ pn++;
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ if(dn == pn) {
|
|
|
|
|
+ reed_solomon_decode(rs
|
|
|
|
|
+ , data_blocks
|
|
|
|
|
+ , block_size
|
|
|
|
|
+ , dec_fec_blocks
|
|
|
|
|
+ , fec_block_nos
|
|
|
|
|
+ , erased_blocks
|
|
|
|
|
+ , dn);
|
|
|
|
|
+ } else {
|
|
|
|
|
+ //error but we continue
|
|
|
|
|
+ err = -1;
|
|
|
|
|
+ }
|
|
|
|
|
+ }
|
|
|
|
|
+ data_blocks += ds;
|
|
|
|
|
+ marks += ds;
|
|
|
|
|
+ fec_blocks += ps;
|
|
|
|
|
+ fec_marks += ps;
|
|
|
|
|
+ }
|
|
|
|
|
+
|
|
|
|
|
+ return err;
|
|
|
|
|
+}
|
|
|
|
|
+
|